Repository: PEtab-dev/petab_select (language: python)

Issue:
Using python >= 3.8 function?
Running tests on python 3.7 gives

```
petab_select/model_space.py:6: in <module>
    from typing import (
E   ImportError: cannot import name 'get_args' from 'typing' (/opt/hostedtoolcache/Python/3.7.12/x64/lib/python3.7/typing.py)
```

(see e.g. https://github.com/PEtab-dev/petab_select/runs/5360619080?check_suite_focus=true)

runs w/o problems on python >= 3.8. So, need to either fix this function use or require >= 3.8 as minimum version.


Discussion / hints:

Fix patch:
diff --git a/.flake8 b/.flake8
new file mode 100644
index 00000000..260dd214
--- /dev/null
+++ b/.flake8
@@ -0,0 +1,9 @@
+########################
+# Flake8 Configuration #
+########################
+
+[flake8]
+
+per-file-ignores =
+    # Imported but unused
+    */__init__.py:F401,F403
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 00000000..6a40e862
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,56 @@
+name: CI
+
+# trigger
+on:
+  push:
+    branches:
+    - main
+  pull_request:
+  schedule:
+    # run Tuesday and Friday at 02:00 UTC
+    - cron: '00 2 * * TUE,FRI'
+
+jobs:
+  base:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        # test on latest and minimum python version
+        python-version: ['3.10', '3.8']
+
+    steps:
+    - name: Check out repository
+      uses: actions/checkout@v2
+
+    - name: Prepare python ${{ matrix.python-version }}
+      uses: actions/setup-python@v1
+      with:
+        python-version: ${{ matrix.python-version }}
+
+    - name: Cache
+      uses: actions/cache@v1
+      with:
+        path: ~ /.cache
+        key: ci-${{ runner.os }}-${{ matrix.python-version }}-base
+
+    - name: Install system dependencies
+      run: |
+        sudo apt-get update
+        sudo apt-get install -y \
+          swig \
+          libatlas-base-dev \
+          libhdf5-serial-dev
+
+    - name: Install Python dependencies
+      run: pip install -r requirements_dev.txt
+
+    - name: Run tox
+      run: python -m tox
+
+    - name: Run pre-commit hooks
+      run: pre-commit run --all-files
+
+    - name: Coverage
+      uses: codecov/codecov-action@v1
+      with:
+        file: ./coverage.xml
diff --git a/.gitignore b/.gitignore
index 596a5ecf..8eb8a830 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,14 @@ __pycache__
 .ipynb_checkpoints
 dist
 amici_models
+*egg-info
+
+# coverage
+.coverage*
+coverage.xml
+
+# tox
+.tox
+
+# osx
+.DS_Store
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 00000000..a3220358
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,46 @@
+# This is run as a precondition to commits, run manually via `pre-commit run`
+
+# When adding new hooks, it may make sense to once run
+# `pre-commit run --all-files` as by default only changed files are checked
+
+repos:
+- repo: https://github.com/psf/black
+  rev: 22.6.0
+  hooks:
+  - id: black
+    description: The uncompromising code formatter
+- repo: https://github.com/pycqa/isort
+  rev: 5.10.1
+  hooks:
+  - id: isort
+    name: isort (python)
+  - id: isort
+    name: isort (cython)
+    types: [cython]
+  - id: isort
+    name: isort (pyi)
+    types: [pyi]
+- repo: https://github.com/nbQA-dev/nbQA
+  rev: 1.3.1
+  hooks:
+    - id: nbqa-black
+    - id: nbqa-pyupgrade
+      args: [--py36-plus]
+    - id: nbqa-isort
+- repo: https://github.com/pre-commit/pre-commit-hooks
+  rev: v4.3.0
+  hooks:
+  - id: check-yaml
+    description: Check yaml files for parseable syntax
+  - id: check-added-large-files
+    description: Prevent large files from being committed
+  - id: check-merge-conflict
+    description: Check for files that contain merge conflict strings
+  - id: check-symlinks
+    description: Check for symlinks which do not point to anything
+  - id: trailing-whitespace
+    description: Trim trailing whitespaces
+  - id: end-of-file-fixer
+    description: Fix empty lines at ends of files
+  - id: detect-private-key
+    description: Detects the presence of private keys
diff --git a/changes.md b/changes.md
index c1f0687e..bbc28d63 100644
--- a/changes.md
+++ b/changes.md
@@ -1,5 +1,60 @@
 # Changes
 
+## 0.1.0
+There are several new candidate spaces, including the FAMoS method. There are also code quality checks (`black` / `isort` / notebook checks) and tests / CI with `tox`.
+
+Below are some significant changes. One major change is that the predecessor model is no longer explicitly specified, instead it is taken as the best model from a set of provided models.
+This means a calibration tool only needs to calibrate all models that PEtab Select provides, then send all calibrated models back to PEtab select, at each iteration.
+This change is implemented to support the FAMoS method.
+
+### `CandidateSpace`
+- new candidate spaces
+  - `BidirectionalCandidateSpace`
+    - tests models with both the forward and backward method at each iteration
+  - `ForwardAndBackwardCandidateSpace`
+    - alternates between the forward and backward method
+  - `LateralCandidateSpace`
+    - moves "sideways" through the model space with parameter swaps
+      - a lateral move is where one parameter is turned off, and another parameter is turned on, simultaneously
+  - `FamosCandidateSpace`
+    - implements the FAMoS method ( https://doi.org/10.1371/journal.pcbi.1007230 )
+- candidate spaces have a `update_after_calibration` method
+  - this provides functionality to support FAMoS operations that occur inbetween calibration iterations
+
+### `Problem`
+- renamed method `add_calibrated_models` to `exclude_models`
+- new method `exclude_model_hashes`
+- cleaned `get_best`
+
+### `ui.candidates`
+- logic change: predecessor model is taken from the following in descending priority, to support the FAMoS method
+  - best of `newly_calibarated_models`
+  - best of `calibrated_models`
+  - `previous_predecessor_model`
+- refactored to instead take newly calibrated models as an argument, instead of assuming they are at `candidate_space.models`
+- renamed `history` to `calibrated_models`
+
+### CLI
+- `candidates`
+  - see logic change in `ui.candidates`
+    - breaking change: renamed `--predecessor` to `--previous-predecessor-model` to ensure this logic change is noticed by developers
+    - removed deprecated `--initial` argument
+    - removed `--best`
+    - added argument `--calibrated-models` to specify file(s) that contain all calibrated models from all iterations so far
+    - added argument `--newly-calibrated-models` to specify file(s) that contain all calibrates models from the latest iteration
+  - renamed arguments to align `cli.candidates` with `ui.candidates`
+    - renamed `--yaml` to `--problem`
+    - renamed `--excluded-model-file` to `--excluded-model-file`
+    - renamed `--excluded-model-hash-file` to `--excluded-model-hashes`
+- `best`
+  - renamed arguments to align `cli.best` with `ui.best`
+    - renamed `--yaml` to `--problem`
+    - renamed `--models_yaml` to `--models`
+- `model_to_petab` and `models_to_petab`, respectively
+  - renamed to align `ui` and `cli` methods
+    - renamed `--yaml` to `model` and `models`, respectively
+    - renamed shorthand `-m` for `--model_id` to `-i` to avoid conflict with new `--model/--models` arguments
+
 ## 0.0.9
 Fixed automatic exclusions from the Python interface. Now, exclusions are only automatically managed in `Problem`s (thereby also the model space and model subspaces) and `CandidateSpace`s.
 
diff --git a/doc/examples/example_cli_famos.ipynb b/doc/examples/example_cli_famos.ipynb
new file mode 100644
index 00000000..945a764c
--- /dev/null
+++ b/doc/examples/example_cli_famos.ipynb
@@ -0,0 +1,209 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "e28947de",
+   "metadata": {},
+   "source": [
+    "# FAMoS in PEtab Select\n",
+    "\n",
+    "In this notebook, the FAMoS algorithm [1] is demonstrated.\n",
+    "\n",
+    "[1] Gabel M, Hohl T, Imle A, Fackler OT, Graw F (2019) FAMoS: A Flexible and dynamic Algorithm for Model Selection to analyse complex systems dynamics. PLOS Computational Biology 15(8): e1007230. https://doi.org/10.1371/journal.pcbi.1007230"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "484b7f3c",
+   "metadata": {},
+   "source": [
+    "The model space contains 65536 models. In normal PEtab Select workflows, a calibration tool would take the place of the `example_cli_famos_calibration_tool.py` script. This script emulates a calibration tool: it takes PEtab Select models and assigns criterion values to them, based on previous calibration results for the same models."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "1f04dce0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Cleanup the state and candidate models output by a previous run of this notebook\n",
+    "import shutil\n",
+    "from pathlib import Path\n",
+    "\n",
+    "from example_cli_famos_helpers import (\n",
+    "    expected_criterion_values,\n",
+    "    parse_summary_to_progress_list,\n",
+    "    petab_select_problem_yaml,\n",
+    ")\n",
+    "\n",
+    "output_path = Path().resolve() / \"output_famos\"\n",
+    "output_path_str = str(output_path)\n",
+    "shutil.rmtree(output_path_str)\n",
+    "output_path.mkdir(exist_ok=False, parents=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "a81560e6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import petab\n",
+    "\n",
+    "import petab_select\n",
+    "from petab_select import ESTIMATE, FamosCandidateSpace, Method, Model\n",
+    "from petab_select.constants import Criterion\n",
+    "from petab_select.model import default_compare\n",
+    "\n",
+    "state = str(output_path / 'state.dill')\n",
+    "\n",
+    "# Each iteration of model selection is described as a 2-tuple here.\n",
+    "# First value is the model selection method.\n",
+    "# Second value are relative change in parameter indices that correspond\n",
+    "# to the best model from this iteration.\n",
+    "# e.g. `(Method.FORWARD, {3})` states that the best model from a forward move\n",
+    "# is the model that now estimates the parameter at index 3.\n",
+    "expected_progress_list = [\n",
+    "    (Method.LATERAL, set()),\n",
+    "    (Method.LATERAL, {4, 15}),\n",
+    "    (Method.LATERAL, {9, 13}),\n",
+    "    (Method.FORWARD, set()),\n",
+    "    (Method.FORWARD, {3}),\n",
+    "    (Method.FORWARD, {11}),\n",
+    "    (Method.BACKWARD, set()),\n",
+    "    (Method.BACKWARD, {6}),\n",
+    "    (Method.BACKWARD, {10}),\n",
+    "    (Method.BACKWARD, {8}),\n",
+    "    (Method.BACKWARD, {14}),\n",
+    "    (Method.BACKWARD, {1}),\n",
+    "    (Method.BACKWARD, {16}),\n",
+    "    (Method.BACKWARD, {4}),\n",
+    "    (Method.FORWARD, set()),\n",
+    "    (Method.LATERAL, set()),\n",
+    "    (Method.MOST_DISTANT, {2, 3, 4, 5, 6, 7, 9, 11, 12, 13, 15}),\n",
+    "    (Method.LATERAL, {16, 7}),\n",
+    "    (Method.LATERAL, {5, 12}),\n",
+    "    (Method.LATERAL, {13, 15}),\n",
+    "    (Method.LATERAL, {1, 6}),\n",
+    "    (Method.FORWARD, set()),\n",
+    "    (Method.FORWARD, {3}),\n",
+    "    (Method.FORWARD, {7}),\n",
+    "    (Method.FORWARD, {2}),\n",
+    "    (Method.FORWARD, {11}),\n",
+    "    (Method.BACKWARD, set()),\n",
+    "    (Method.BACKWARD, {7}),\n",
+    "    (Method.BACKWARD, {16}),\n",
+    "    (Method.BACKWARD, {4}),\n",
+    "    (Method.FORWARD, set()),\n",
+    "    (Method.LATERAL, set()),\n",
+    "    (Method.LATERAL, {9, 15}),\n",
+    "    (Method.FORWARD, set()),\n",
+    "    (Method.BACKWARD, set()),\n",
+    "    (Method.LATERAL, set()),\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7202f6c6",
+   "metadata": {},
+   "source": [
+    "The predecessor model is some model from the model space, and is defined in the PEtab Select problem YAML file."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "bb1a5144",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%bash -s \"$petab_select_problem_yaml\" \"$output_path_str\"\n",
+    "\n",
+    "petab_select_problem_yaml=$1\n",
+    "output_path_str=$2\n",
+    "\n",
+    "problem=$petab_select_problem_yaml\n",
+    "state=$output_path_str/state.dill\n",
+    "\n",
+    "newly_calibrated_models=\"\"\n",
+    "calibrated_models=\"\"\n",
+    "empty_output=$(echo -e \"null\\n...\")\n",
+    "\n",
+    "for i in {1..40}\n",
+    "do\n",
+    "\n",
+    "output=$output_path_str/models_$i.yaml\n",
+    "calibrated_output=$output_path_str/calibrated_models_$i.yaml\n",
+    "\n",
+    "petab_select candidates \\\n",
+    "--problem $problem \\\n",
+    "--state $state \\\n",
+    "--output $output \\\n",
+    "$newly_calibrated_models \\\n",
+    "$calibrated_models \\\n",
+    "--relative-paths\n",
+    "\n",
+    "# Replace this line with a tool that calibrates the models.\n",
+    "# The script also changes model IDs for easier analysis in this example.\n",
+    "python example_cli_famos_calibration_tool.py $output $calibrated_output\n",
+    "\n",
+    "if [ \"$(cat $calibrated_output)\" == \"$empty_output\" ]\n",
+    "then\n",
+    "        # End the model selection if no models were provided by PEtab Select.\n",
+    "        break\n",
+    "else\n",
+    "        newly_calibrated_models=\"--newly-calibrated-models \"$calibrated_output\n",
+    "        calibrated_models+=\"--calibrated-models \"$calibrated_output\" \"\n",
+    "fi\n",
+    "\n",
+    "done"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "93caf071",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "progress_list = parse_summary_to_progress_list('output_famos/summary.tsv')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "cb61d0f7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert progress_list == expected_progress_list"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/doc/examples/example_cli_famos_calibration_tool.py b/doc/examples/example_cli_famos_calibration_tool.py
new file mode 100644
index 00000000..1f589fba
--- /dev/null
+++ b/doc/examples/example_cli_famos_calibration_tool.py
@@ -0,0 +1,32 @@
+import sys
+
+import pandas as pd
+from example_cli_famos_helpers import calibrate, set_model_id
+
+import petab_select
+from petab_select import ESTIMATE, Criterion, Model
+
+models_yaml = sys.argv[1]
+calibrated_models_yaml = sys.argv[2]
+
+models = petab_select.model.models_from_yaml_list(models_yaml)
+predecessor_model_hashes = set()
+for model in models:
+    set_model_id(model=model)
+    calibrate(model=model)
+    predecessor_model_hashes |= {model.predecessor_model_hash}
+petab_select.model.models_to_yaml_list(
+    models=models, output_yaml=calibrated_models_yaml
+)
+
+model_ids = sorted(model.model_id for model in models)
+
+if len(predecessor_model_hashes) == 0:
+    pass
+elif len(predecessor_model_hashes) == 1:
+    (predecessor_model_hash,) = predecessor_model_hashes
+else:
+    print(
+        'The models of this iteration somehow have different predecessor models.\n'
+        + '\n'.join(predecessor_model_hashes)
+    )
diff --git a/doc/examples/example_cli_famos_helpers.py b/doc/examples/example_cli_famos_helpers.py
new file mode 100644
index 00000000..cdbd2cdd
--- /dev/null
+++ b/doc/examples/example_cli_famos_helpers.py
@@ -0,0 +1,77 @@
+from pathlib import Path
+from typing import Tuple
+
+import pandas as pd
+
+import petab_select
+from petab_select import ESTIMATE, Method, Model
+
+input_path = (
+    Path(__file__).resolve().parent.parent.parent
+    / "test"
+    / "candidate_space"
+    / "input"
+    / "famos_synthetic"
+)
+petab_select_problem_yaml = str(
+    input_path / "select" / "FAMoS_2019_petab_select_problem.yaml"
+)
+expected_criterion_values = dict(
+    pd.read_csv(
+        input_path / "test_files" / "calibration_results.tsv", sep="\t"
+    ).set_index('model_id')['AICc']
+)
+
+petab_select_problem = petab_select.Problem.from_yaml(
+    petab_select_problem_yaml
+)
+criterion = petab_select_problem.criterion
+
+
+def set_model_id(model: Model) -> None:
+    """Set the model ID to a binary string for easier analysis."""
+    model.model_id = "M_" + ''.join(
+        str(v) for v in model.model_subspace_indices
+    )
+
+
+def calibrate(
+    model: Model,
+    criterion=criterion,
+    expected_criterion_values=expected_criterion_values,
+) -> None:
+    """Set the criterion value for a model."""
+    model.set_criterion(
+        criterion=criterion,
+        value=float(expected_criterion_values[model.model_id]),
+    )
+
+
+def parse_summary_to_progress_list(summary_tsv: str) -> Tuple[Method, set]:
+    """Get progress information from the summary file."""
+    df_raw = pd.read_csv(summary_tsv, sep='\t')
+    df = df_raw.loc[~pd.isnull(df_raw["predecessor change"])]
+
+    parameter_list = list(
+        petab_select_problem.model_space.model_subspaces[
+            'model_subspace_1'
+        ].parameters
+    )
+
+    progress_list = []
+
+    for index, (_, row) in enumerate(df.iterrows()):
+        method = Method(row["method"])
+
+        model = {
+            1 + parameter_list.index(parameter_id)
+            for parameter_id in eval(row["current model"])
+        }
+        if index == 0:
+            model0 = model
+
+        difference = model.symmetric_difference(model0)
+        progress_list.append((method, difference))
+        model0 = model
+
+    return progress_list
diff --git a/doc/examples/model_selection/calibrated/M1_4_0.yaml b/doc/examples/model_selection/calibrated_M1_4.yaml
similarity index 72%
rename from doc/examples/model_selection/calibrated/M1_4_0.yaml
rename to doc/examples/model_selection/calibrated_M1_4.yaml
index e96df2a2..90de1f16 100644
--- a/doc/examples/model_selection/calibrated/M1_4_0.yaml
+++ b/doc/examples/model_selection/calibrated_M1_4.yaml
@@ -8,4 +8,4 @@ parameters:
 estimated_parameters:
   k2: 0.15
   k3: 0.0
-petab_yaml: ../petab_problem.yaml
+petab_yaml: ../model_selection/petab_problem.yaml
diff --git a/doc/examples/model_selection/calibrated/M1_7_0.yaml b/doc/examples/model_selection/calibrated_M1_7.yaml
similarity index 74%
rename from doc/examples/model_selection/calibrated/M1_7_0.yaml
rename to doc/examples/model_selection/calibrated_M1_7.yaml
index ac01f6a8..fdcfc217 100644
--- a/doc/examples/model_selection/calibrated/M1_7_0.yaml
+++ b/doc/examples/model_selection/calibrated_M1_7.yaml
@@ -9,4 +9,4 @@ estimated_parameters:
   k1: 0.25
   k2: 0.1
   k3: 0.0
-petab_yaml: ../petab_problem.yaml
+petab_yaml: ../model_selection/petab_problem.yaml
diff --git a/doc/examples/model_selection/calibrated/calibrated_first_iteration.yaml b/doc/examples/model_selection/calibrated_first_iteration.yaml
similarity index 61%
rename from doc/examples/model_selection/calibrated/calibrated_first_iteration.yaml
rename to doc/examples/model_selection/calibrated_first_iteration.yaml
index 7e4c66e9..97d024de 100644
--- a/doc/examples/model_selection/calibrated/calibrated_first_iteration.yaml
+++ b/doc/examples/model_selection/calibrated_first_iteration.yaml
@@ -1,31 +1,29 @@
 - criteria:
     AIC: 180
-  model_id: M1_0_0
+  model_id: M1_0
   parameters:
     k1: 0.0
     k2: 0.0
     k3: 0.0
   estimated_parameters: null
-  petab_yaml: ../petab_problem.yaml
+  petab_yaml: ../model_selection/petab_problem.yaml
 - criteria:
     AIC: 100
-  model_id: M1_1_0
+  model_id: M1_1
   parameters:
     k1: 0.2
     k2: 0.1
     k3: estimate
   estimated_parameters:
     k3: 0.0
-  petab_yaml: ../petab_problem.yaml
+  petab_yaml: ../model_selection/petab_problem.yaml
 - criteria:
     AIC: 50
-  model_id: M1_2_0
+  model_id: M1_2
   parameters:
     k1: 0.2
     k2: estimate
     k3: 0.0
   estimated_parameters:
     k2: 0.05
-  petab_yaml: ../petab_problem.yaml
-
-
+  petab_yaml: ../model_selection/petab_problem.yaml
diff --git a/doc/examples/model_selection/observables.tsv b/doc/examples/model_selection/observables.tsv
index ab25269a..b244d788 100644
--- a/doc/examples/model_selection/observables.tsv
+++ b/doc/examples/model_selection/observables.tsv
@@ -1,2 +1,2 @@
 observableId	observableFormula	noiseFormula
-obs_x2	x2+k4	noiseParameter1_obs_x2
+obs_x2	x2+k4+k5	noiseParameter1_obs_x2
diff --git a/doc/examples/model_selection/parameters.tsv b/doc/examples/model_selection/parameters.tsv
index 7a25c2ed..04850368 100644
--- a/doc/examples/model_selection/parameters.tsv
+++ b/doc/examples/model_selection/parameters.tsv
@@ -3,4 +3,5 @@ k1	k_{1}	lin	0	1e3	0.2	1
 k2	k_{2}	lin	0	1e3	0.1	1
 k3	k_{3}	lin	0	1e3	0	1
 k4	k_{4}	lin	0	1e3	0	0
+k5	k_{5}	lin	0	1e3	0	0
 sigma_x2	\sigma_{x2}	lin	1e-5	1e3	0.15	1
diff --git a/doc/examples/workflow_cli.ipynb b/doc/examples/workflow_cli.ipynb
index 5067a2cd..a2bc1c1c 100644
--- a/doc/examples/workflow_cli.ipynb
+++ b/doc/examples/workflow_cli.ipynb
@@ -67,24 +67,15 @@
    "execution_count": 2,
    "id": "eab391ee",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/home/dilan/Documents/future_annex/model_selection/packages/petab_select/petab_select/model_subspace.py:236: UserWarning: The candidate space has already reached its limit of accepted models.\n",
-      "  warnings.warn('The candidate space has already reached its limit of accepted models.')  # noqa: E501\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "%%bash\n",
     "\n",
     "petab_select candidates \\\n",
-    "-y model_selection/petab_select_problem.yaml \\\n",
-    "-s output/state.dill \\\n",
-    "-o output/models_1.yaml \\\n",
-    "-m brute_force \\\n",
+    "--problem model_selection/petab_select_problem.yaml \\\n",
+    "--state output/state.dill \\\n",
+    "--output output/models_1.yaml \\\n",
+    "--method brute_force \\\n",
     "--relative-paths \\\n",
     "-l 3"
    ]
@@ -157,7 +148,7 @@
     }
    ],
    "source": [
-    "with open('output/models_1.yaml', 'r') as f:\n",
+    "with open('output/models_1.yaml') as f:\n",
     "    print(f.read())"
    ]
   },
@@ -174,7 +165,7 @@
    "id": "d8008447",
    "metadata": {},
    "source": [
-    "Between iterations, the models from the first iteration have been calibrated, and the model with the best criterion value is `M1_2_0`. Here, PEtab Select will identify `M1_2_0` as the best model from a YAML file of calibrated models. In the following iterations, a specific model will be provided."
+    "Between iterations, the models from the first iteration have been calibrated, and the model with the best criterion value is `M1_2`. Here, PEtab Select will identify `M1_2` as the best model from a YAML file of calibrated models. In the following iterations, a specific model will be provided."
    ]
   },
   {
@@ -189,41 +180,39 @@
      "text": [
       "- criteria:\n",
       "    AIC: 180\n",
-      "  model_id: M1_0_0\n",
+      "  model_id: M1_0\n",
       "  parameters:\n",
       "    k1: 0.0\n",
       "    k2: 0.0\n",
       "    k3: 0.0\n",
       "  estimated_parameters: null\n",
-      "  petab_yaml: ../petab_problem.yaml\n",
+      "  petab_yaml: ../model_selection/petab_problem.yaml\n",
       "- criteria:\n",
       "    AIC: 100\n",
-      "  model_id: M1_1_0\n",
+      "  model_id: M1_1\n",
       "  parameters:\n",
       "    k1: 0.2\n",
       "    k2: 0.1\n",
       "    k3: estimate\n",
       "  estimated_parameters:\n",
       "    k3: 0.0\n",
-      "  petab_yaml: ../petab_problem.yaml\n",
+      "  petab_yaml: ../model_selection/petab_problem.yaml\n",
       "- criteria:\n",
       "    AIC: 50\n",
-      "  model_id: M1_2_0\n",
+      "  model_id: M1_2\n",
       "  parameters:\n",
       "    k1: 0.2\n",
       "    k2: estimate\n",
       "    k3: 0.0\n",
       "  estimated_parameters:\n",
       "    k2: 0.05\n",
-      "  petab_yaml: ../petab_problem.yaml\n",
-      "\n",
-      "\n",
+      "  petab_yaml: ../model_selection/petab_problem.yaml\n",
       "\n"
      ]
     }
    ],
    "source": [
-    "with open('model_selection/calibrated/calibrated_first_iteration.yaml', 'r') as f:\n",
+    "with open('model_selection/calibrated_first_iteration.yaml') as f:\n",
     "    print(f.read())"
    ]
   },
@@ -245,11 +234,11 @@
     "%%bash\n",
     "\n",
     "petab_select candidates \\\n",
-    "-y model_selection/petab_select_problem.yaml \\\n",
-    "-s output/state.dill \\\n",
-    "-o output/models_2.yaml \\\n",
-    "-m forward \\\n",
-    "-b model_selection/calibrated/calibrated_first_iteration.yaml \\\n",
+    "--problem model_selection/petab_select_problem.yaml \\\n",
+    "--state output/state.dill \\\n",
+    "--output output/models_2.yaml \\\n",
+    "--method forward \\\n",
+    "--newly-calibrated-models model_selection/calibrated_first_iteration.yaml \\\n",
     "--relative-paths"
    ]
   },
@@ -258,7 +247,7 @@
    "id": "ac1f8a5f",
    "metadata": {},
    "source": [
-    "`M1_2_0` has one estimated parameter, `k2` [1]. As expected, the new candidates identified with the forward method have two estimated parameters, and one of them is `k2`.\n",
+    "`M1_2` has one estimated parameter, `k2` [1]. As expected, the new candidates identified with the forward method have two estimated parameters, and one of them is `k2`.\n",
     "\n",
     "[1] There may be additional estimated parameters specified in the PEtab problem."
    ]
@@ -287,7 +276,7 @@
       "    k2: estimate\n",
       "    k3: estimate\n",
       "  petab_yaml: ../model_selection/petab_problem.yaml\n",
-      "  predecessor_model_hash: b91dd6d1144253c4462597de1fff07cb20c6e45ab33d55a8282bf2b3df25d1226ad19baf488a96852a626df1b56e1e618929f0c0627ed6f2c3f69ebb0bcbb6fc\n",
+      "  predecessor_model_hash: c95cc8f9f2b945c996856212e7d49784f72cc429dbdbb9bc04adad610f1a8103a58298c714d0518f16a5a002db5a9c27847f5da6a863637178a15a15f2fc294d\n",
       "- criteria: {}\n",
       "  estimated_parameters: {}\n",
       "  model_hash: db8700c079c8347123adc89b7f5112256c4aaebd2af0f6e32e7582f398b2c1e5e85e588cdcc56bab054c001b96a9b42b02174266927f879d7f78e8ac5d2c33e6\n",
@@ -302,13 +291,13 @@
       "    k2: estimate\n",
       "    k3: 0\n",
       "  petab_yaml: ../model_selection/petab_problem.yaml\n",
-      "  predecessor_model_hash: b91dd6d1144253c4462597de1fff07cb20c6e45ab33d55a8282bf2b3df25d1226ad19baf488a96852a626df1b56e1e618929f0c0627ed6f2c3f69ebb0bcbb6fc\n",
+      "  predecessor_model_hash: c95cc8f9f2b945c996856212e7d49784f72cc429dbdbb9bc04adad610f1a8103a58298c714d0518f16a5a002db5a9c27847f5da6a863637178a15a15f2fc294d\n",
       "\n"
      ]
     }
    ],
    "source": [
-    "with open('output/models_2.yaml', 'r') as f:\n",
+    "with open('output/models_2.yaml') as f:\n",
     "    print(f.read())"
    ]
   },
@@ -325,7 +314,7 @@
    "id": "afba880c",
    "metadata": {},
    "source": [
-    "Models `M1_4_0` and `M1_6_0` were calibrated, with `M1_4_0` having the superior criterion value, so is used as the predecessor model in the next forward search."
+    "Models `M1_4` and `M1_6` were calibrated, with `M1_4` having the superior criterion value, so is used as the predecessor model in the next forward search."
    ]
   },
   {
@@ -348,13 +337,13 @@
       "estimated_parameters:\n",
       "  k2: 0.15\n",
       "  k3: 0.0\n",
-      "petab_yaml: ../petab_problem.yaml\n",
+      "petab_yaml: ../model_selection/petab_problem.yaml\n",
       "\n"
      ]
     }
    ],
    "source": [
-    "with open('model_selection/calibrated/M1_4_0.yaml', 'r') as f:\n",
+    "with open('model_selection/calibrated_M1_4.yaml') as f:\n",
     "    print(f.read())"
    ]
   },
@@ -368,11 +357,11 @@
     "%%bash\n",
     "\n",
     "petab_select candidates \\\n",
-    "-y model_selection/petab_select_problem.yaml \\\n",
-    "-s output/state.dill \\\n",
-    "-o output/models_3.yaml \\\n",
-    "-m forward \\\n",
-    "-p model_selection/calibrated/M1_4_0.yaml \\\n",
+    "--problem model_selection/petab_select_problem.yaml \\\n",
+    "--state output/state.dill \\\n",
+    "--output output/models_3.yaml \\\n",
+    "--method forward \\\n",
+    "--newly-calibrated-models model_selection/calibrated_M1_4.yaml \\\n",
     "--relative-paths"
    ]
   },
@@ -408,13 +397,13 @@
       "    k2: estimate\n",
       "    k3: estimate\n",
       "  petab_yaml: ../model_selection/petab_problem.yaml\n",
-      "  predecessor_model_hash: 3a5dadef211388a3da6610b1c057d16a135dcddd5273efa88b1875bbafb153afeb4398a8ed02cf04e0fa7c04348265bebe2b02b8c0def80e2fb6f329e469e7f7\n",
+      "  predecessor_model_hash: 0e0f00ee37d03db0f6ad7bb9cb92eb4ae3aef85fa1e1d9d38dfc5954fbcf62e8768803bedd42fb3638b681d4840e352e5bfc29a1800247ebeaddfbdb4baca46b\n",
       "\n"
      ]
     }
    ],
    "source": [
-    "with open('output/models_3.yaml', 'r') as f:\n",
+    "with open('output/models_3.yaml') as f:\n",
     "    print(f.read())"
    ]
   },
@@ -448,13 +437,13 @@
       "  k1: 0.25\n",
       "  k2: 0.1\n",
       "  k3: 0.0\n",
-      "petab_yaml: ../petab_problem.yaml\n",
+      "petab_yaml: ../model_selection/petab_problem.yaml\n",
       "\n"
      ]
     }
    ],
    "source": [
-    "with open('model_selection/calibrated/M1_7_0.yaml', 'r') as f:\n",
+    "with open('model_selection/calibrated_M1_7.yaml') as f:\n",
     "    print(f.read())"
    ]
   },
@@ -468,11 +457,11 @@
     "%%bash\n",
     "\n",
     "petab_select candidates \\\n",
-    "-y model_selection/petab_select_problem.yaml \\\n",
-    "-s output/state.dill \\\n",
-    "-o output/models_4.yaml \\\n",
-    "-m forward \\\n",
-    "-p model_selection/calibrated/M1_7_0.yaml \\\n",
+    "--problem model_selection/petab_select_problem.yaml \\\n",
+    "--state output/state.dill \\\n",
+    "--output output/models_4.yaml \\\n",
+    "--method forward \\\n",
+    "--newly-calibrated-models model_selection/calibrated_M1_7.yaml \\\n",
     "--relative-paths"
    ]
   },
@@ -493,7 +482,7 @@
     }
    ],
    "source": [
-    "with open('output/models_4.yaml', 'r') as f:\n",
+    "with open('output/models_4.yaml') as f:\n",
     "    print(f.read())"
    ]
   },
@@ -503,7 +492,7 @@
    "metadata": {},
    "source": [
     "## Fifth iteration\n",
-    "Although no additional models are found with a forward search initialized at the best model so far (`M1_7_0`), there are additional models in the model space that are yet to be calibrated, which can be identified with the brute force method."
+    "Although no additional models are found with a forward search initialized at the best model so far (`M1_7`), there are additional models in the model space that are yet to be calibrated, which can be identified with the brute force method."
    ]
   },
   {
@@ -516,14 +505,14 @@
     "%%bash\n",
     "\n",
     "petab_select candidates \\\n",
-    "-y model_selection/petab_select_problem.yaml \\\n",
-    "-s output/state.dill \\\n",
-    "-o output/models.yaml \\\n",
-    "-m brute_force \\\n",
-    "-e output/models_1.yaml \\\n",
-    "-e output/models_2.yaml \\\n",
-    "-e output/models_3.yaml \\\n",
-    "-e output/models_4.yaml \\\n",
+    "--problem model_selection/petab_select_problem.yaml \\\n",
+    "--state output/state.dill \\\n",
+    "--output output/models.yaml \\\n",
+    "--method brute_force \\\n",
+    "--excluded-models output/models_1.yaml \\\n",
+    "--excluded-models output/models_2.yaml \\\n",
+    "--excluded-models output/models_3.yaml \\\n",
+    "--excluded-models output/models_4.yaml \\\n",
     "--relative-paths"
    ]
   },
@@ -572,7 +561,7 @@
     }
    ],
    "source": [
-    "with open('output/models.yaml', 'r') as f:\n",
+    "with open('output/models.yaml') as f:\n",
     "    print(f.read())"
    ]
   },
@@ -595,10 +584,12 @@
     "%%bash\n",
     "\n",
     "petab_select best \\\n",
-    "-y model_selection/petab_select_problem.yaml \\\n",
-    "-m model_selection/calibrated/calibrated_first_iteration.yaml \\\n",
-    "-o output/best_model.yaml \\\n",
-    "-s output/state.dill \\\n",
+    "--problem model_selection/petab_select_problem.yaml \\\n",
+    "--models model_selection/calibrated_first_iteration.yaml \\\n",
+    "--models model_selection/calibrated_M1_4.yaml \\\n",
+    "--models model_selection/calibrated_M1_7.yaml \\\n",
+    "--output output/best_model.yaml \\\n",
+    "--state output/state.dill \\\n",
     "--relative-paths"
    ]
   },
@@ -613,17 +604,19 @@
      "output_type": "stream",
      "text": [
       "criteria:\n",
-      "  AIC: 50\n",
+      "  AIC: 20\n",
       "estimated_parameters:\n",
-      "  k2: 0.05\n",
+      "  k1: 0.25\n",
+      "  k2: 0.1\n",
+      "  k3: 0.0\n",
       "model_hash: null\n",
-      "model_id: M1_2_0\n",
+      "model_id: M1_7_0\n",
       "model_subspace_id: null\n",
       "model_subspace_indices: null\n",
       "parameters:\n",
-      "  k1: 0.2\n",
+      "  k1: estimate\n",
       "  k2: estimate\n",
-      "  k3: 0.0\n",
+      "  k3: estimate\n",
       "petab_yaml: ../model_selection/petab_problem.yaml\n",
       "predecessor_model_hash: null\n",
       "\n"
@@ -631,7 +624,7 @@
     }
    ],
    "source": [
-    "with open('output/best_model.yaml', 'r') as f:\n",
+    "with open('output/best_model.yaml') as f:\n",
     "    print(f.read())"
    ]
   },
@@ -661,8 +654,8 @@
     "%%bash\n",
     "\n",
     "petab_select model_to_petab \\\n",
-    "-y output/best_model.yaml \\\n",
-    "-o output/best_model_petab"
+    "--model output/best_model.yaml \\\n",
+    "--output output/best_model_petab"
    ]
   }
  ],
diff --git a/doc/examples/workflow_python.ipynb b/doc/examples/workflow_python.ipynb
index 0f2feaa5..361f12a6 100644
--- a/doc/examples/workflow_python.ipynb
+++ b/doc/examples/workflow_python.ipynb
@@ -46,14 +46,12 @@
    ],
    "source": [
     "import petab_select\n",
-    "from petab_select import (\n",
-    "    Model,\n",
-    "    ForwardCandidateSpace,\n",
-    ")\n",
-    "\n",
+    "from petab_select import ForwardCandidateSpace, Model\n",
     "\n",
     "# Load the PEtab Select problem.\n",
-    "select_problem = petab_select.Problem.from_yaml('model_selection/petab_select_problem.yaml')\n",
+    "select_problem = petab_select.Problem.from_yaml(\n",
+    "    'model_selection/petab_select_problem.yaml'\n",
+    ")\n",
     "# Fake criterion values as a surrogate for a model calibration tool.\n",
     "fake_criterion = {\n",
     "    'M1_0': 200,\n",
@@ -70,7 +68,7 @@
     "def print_model(model: Model) -> None:\n",
     "    \"\"\"Helper method to view model attributes.\"\"\"\n",
     "    print(\n",
-    "f\"\"\"\\\n",
+    "        f\"\"\"\\\n",
     "Model subspace ID: {model.model_subspace_id}\n",
     "PEtab YAML location: {model.petab_yaml}\n",
     "Custom model parameters: {model.parameters}\n",
@@ -83,15 +81,21 @@
     "\n",
     "def calibrate(model: Model, fake_criterion=fake_criterion) -> None:\n",
     "    \"\"\"Set model criterion values to fake values that could be the output of a calibration tool.\n",
-    "    \n",
+    "\n",
     "    Each model subspace in this problem contains only one model, so a model-specific criterion can\n",
     "    be indexed by the model subspace ID.\n",
     "    \"\"\"\n",
-    "    model.set_criterion(select_problem.criterion, fake_criterion[model.model_subspace_id])\n",
+    "    model.set_criterion(\n",
+    "        select_problem.criterion, fake_criterion[model.model_subspace_id]\n",
+    "    )\n",
     "\n",
     "\n",
+    "# Calibrated and newly calibrated models should be tracked between iterations.\n",
+    "calibrated_models = {}\n",
+    "newly_calibrated_models = {}\n",
+    "\n",
     "print(\n",
-    "f\"\"\"Information about the model selection problem.\n",
+    "    f\"\"\"Information about the model selection problem.\n",
     "\n",
     "YAML path: {select_problem.yaml_path}\n",
     "Method: {select_problem.method}\n",
@@ -198,7 +202,13 @@
     "# Set fake criterion values that might be the output of a model calibration tool.\n",
     "for candidate_model in candidate_space.models:\n",
     "    calibrate(candidate_model)\n",
-    "select_problem.add_calibrated_models(candidate_space.models)"
+    "\n",
+    "newly_calibrated_models = {\n",
+    "    model.get_hash(): model for model in candidate_space.models\n",
+    "}\n",
+    "calibrated_models.update(newly_calibrated_models)\n",
+    "\n",
+    "select_problem.exclude_models(newly_calibrated_models.values())"
    ]
   },
   {
@@ -222,7 +232,7 @@
     }
    ],
    "source": [
-    "local_best_model = select_problem.get_best(candidate_space.models)\n",
+    "local_best_model = select_problem.get_best(newly_calibrated_models.values())\n",
     "print_model(local_best_model)"
    ]
   },
@@ -253,7 +263,7 @@
     "petab_select.ui.candidates(\n",
     "    problem=select_problem,\n",
     "    candidate_space=candidate_space,\n",
-    "    predecessor_model=select_problem.get_best(candidate_space.models),\n",
+    "    newly_calibrated_models=newly_calibrated_models,\n",
     ");"
    ]
   },
@@ -306,7 +316,13 @@
     "# Set fake criterion values that might be the output of a model calibration tool.\n",
     "for candidate_model in candidate_space.models:\n",
     "    calibrate(candidate_model)\n",
-    "select_problem.add_calibrated_models(candidate_space.models)"
+    "\n",
+    "newly_calibrated_models = {\n",
+    "    model.get_hash(): model for model in candidate_space.models\n",
+    "}\n",
+    "calibrated_models.update(newly_calibrated_models)\n",
+    "\n",
+    "select_problem.exclude_models(newly_calibrated_models.values())"
    ]
   },
   {
@@ -330,7 +346,7 @@
     }
    ],
    "source": [
-    "local_best_model = select_problem.get_best(candidate_space.models)\n",
+    "local_best_model = select_problem.get_best(newly_calibrated_models.values())\n",
     "print_model(local_best_model)"
    ]
   },
@@ -352,7 +368,7 @@
     "petab_select.ui.candidates(\n",
     "    problem=select_problem,\n",
     "    candidate_space=candidate_space,\n",
-    "    predecessor_model=select_problem.get_best(candidate_space.models),\n",
+    "    newly_calibrated_models=newly_calibrated_models,\n",
     ");"
    ]
   },
@@ -398,7 +414,13 @@
     "# Set fake criterion values that might be the output of a model calibration tool.\n",
     "for candidate_model in candidate_space.models:\n",
     "    calibrate(candidate_model)\n",
-    "select_problem.add_calibrated_models(candidate_space.models)"
+    "\n",
+    "newly_calibrated_models = {\n",
+    "    model.get_hash(): model for model in candidate_space.models\n",
+    "}\n",
+    "calibrated_models.update(newly_calibrated_models)\n",
+    "\n",
+    "select_problem.exclude_models(newly_calibrated_models.values())"
    ]
   },
   {
@@ -422,7 +444,7 @@
     }
    ],
    "source": [
-    "local_best_model = select_problem.get_best(candidate_space.models)\n",
+    "local_best_model = select_problem.get_best(newly_calibrated_models.values())\n",
     "print_model(local_best_model)"
    ]
   },
@@ -444,7 +466,7 @@
     "petab_select.ui.candidates(\n",
     "    problem=select_problem,\n",
     "    candidate_space=candidate_space,\n",
-    "    predecessor_model=select_problem.get_best(candidate_space.models),\n",
+    "    newly_calibrated_models=newly_calibrated_models,\n",
     ");"
    ]
   },
@@ -483,7 +505,13 @@
     "# Set fake criterion values that might be the output of a model calibration tool.\n",
     "for candidate_model in candidate_space.models:\n",
     "    calibrate(candidate_model)\n",
-    "select_problem.add_calibrated_models(candidate_space.models)"
+    "\n",
+    "newly_calibrated_models = {\n",
+    "    model.get_hash(): model for model in candidate_space.models\n",
+    "}\n",
+    "calibrated_models.update(newly_calibrated_models)\n",
+    "\n",
+    "select_problem.exclude_models(newly_calibrated_models.values())"
    ]
   },
   {
@@ -507,7 +535,7 @@
     }
    ],
    "source": [
-    "local_best_model = select_problem.get_best(candidate_space.models)\n",
+    "local_best_model = select_problem.get_best(newly_calibrated_models.values())\n",
     "print_model(local_best_model)"
    ]
   },
@@ -529,7 +557,6 @@
     "petab_select.ui.candidates(\n",
     "    problem=select_problem,\n",
     "    candidate_space=candidate_space,\n",
-    "    predecessor_model=select_problem.get_best(candidate_space.models),\n",
     ");"
    ]
   },
@@ -551,7 +578,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Number of candidate models: 0.\n"
+      "Number of candidate models: 1.\n"
      ]
     }
    ],
@@ -588,7 +615,7 @@
     }
    ],
    "source": [
-    "best_model = select_problem.get_best(select_problem.calibrated_models)\n",
+    "best_model = select_problem.get_best(calibrated_models.values())\n",
     "print_model(best_model)"
    ]
   },
@@ -612,7 +639,7 @@
     "petab_select.ui.candidates(\n",
     "    problem=select_problem,\n",
     "    candidate_space=candidate_space,\n",
-    "    excluded_models=select_problem.calibrated_models,\n",
+    "    calibrated_models=calibrated_models,\n",
     ");"
    ]
   },
diff --git a/petab_select/__init__.py b/petab_select/__init__.py
index 20c60b04..bd2994f5 100644
--- a/petab_select/__init__.py
+++ b/petab_select/__init__.py
@@ -1,11 +1,11 @@
 """Model selection extension for PEtab."""
 
-from .candidate_space import *  # noqa: F401, F403
-from .constants import *  # noqa: F401, F403
-from .criteria import *  # noqa: F401, F403
-from .misc import *  # noqa: F401, F403
-from .model import *  # noqa: F401, F403
-from .model_space import *  # noqa: F401, F403
-from .model_subspace import *  # noqa: F401, F403
-from .problem import *  # noqa: F401, F403
-from .ui import *  # noqa: F401, F403
+from .candidate_space import *
+from .constants import *
+from .criteria import *
+from .misc import *
+from .model import *
+from .model_space import *
+from .model_subspace import *
+from .problem import *
+from .ui import *
diff --git a/petab_select/candidate_space.py b/petab_select/candidate_space.py
index 9175b2ae..8fd4b3d8 100644
--- a/petab_select/candidate_space.py
+++ b/petab_select/candidate_space.py
@@ -1,16 +1,34 @@
 """Classes and methods related to candidate spaces."""
 import abc
+import bisect
+import copy
+import csv
+import logging
+import os.path
 import warnings
-from typing import Any, Dict, List, Optional, Union
+from pathlib import Path
+from typing import Any, Callable, Dict, List, Optional, Union
 
 import numpy as np
 from more_itertools import one
 
-from .constants import (ESTIMATE, VIRTUAL_INITIAL_MODEL,
-                        VIRTUAL_INITIAL_MODEL_METHODS, Method)
+from .constants import (
+    ESTIMATE,
+    METHOD,
+    METHOD_SCHEME,
+    MODELS,
+    NEXT_METHOD,
+    PREDECESSOR_MODEL,
+    PREVIOUS_METHODS,
+    TYPE_PATH,
+    VIRTUAL_INITIAL_MODEL,
+    VIRTUAL_INITIAL_MODEL_METHODS,
+    Criterion,
+    Method,
+)
 from .handlers import TYPE_LIMIT, LimitHandler
 from .misc import snake_case_to_camel_case
-from .model import Model
+from .model import Model, default_compare
 
 
 class CandidateSpace(abc.ABC):
@@ -34,25 +52,102 @@ class CandidateSpace(abc.ABC):
             are added to `exclusions`.
         limit:
             A handler to limit the number of accepted models.
+        method:
+            The model selection method of the candidate space.
+        governing_method:
+            Used to store the search method that governs the choice of method during
+            a search. In some cases, this is always the same as the method attribute.
+            An example of a difference is in the bidirectional method, where `governing_method`
+            stores the bidirectional method, whereas `method` may also store the forward or
+            backward methods.
+        retry_model_space_search_if_no_models:
+            Whether a search with a candidate space should be repeated upon failure.
+            Useful for the `BidirectionalCandidateSpace`, which switches directions
+            upon failure.
+        summary_tsv:
+            A string or `pathlib.Path`. A summary of the model selection progress
+            will be written to this file.
         #limited:
         #    A descriptor that handles the limit on the number of accepted models.
         #limit:
         #    Models will fail `self.consider` if `len(self.models) >= limit`.
     """
 
+    governing_method: Method = None
+    method: Method = None
+    retry_model_space_search_if_no_models: bool = False
+
     def __init__(
         self,
         # TODO add MODEL_TYPE = Union[str, Model], str for VIRTUAL_INITIAL_MODEL
         predecessor_model: Optional[Model] = None,
         exclusions: Optional[List[Any]] = None,
         limit: TYPE_LIMIT = np.inf,
+        summary_tsv: TYPE_PATH = None,
     ):
         self.limit = LimitHandler(
             current=self.n_accepted,
             limit=limit,
         )
+        # Each candidate class specifies this as a class attribute.
+        if self.governing_method is None:
+            self.governing_method = self.method
         self.reset(predecessor_model=predecessor_model, exclusions=exclusions)
 
+        self.summary_tsv = summary_tsv
+        if self.summary_tsv is not None:
+            self.summary_tsv = Path(self.summary_tsv)
+            self._setup_summary_tsv()
+
+    def write_summary_tsv(self, row):
+        if self.summary_tsv is None:
+            return
+
+        # Format single values to be valid rows
+        if not isinstance(row, list):
+            row = [
+                row,
+                *([''] * 5),
+            ]
+
+        with open(self.summary_tsv, 'a', encoding="utf-8") as f:
+            writer = csv.writer(f, delimiter='\t')
+            writer.writerow(row)
+
+    def _setup_summary_tsv(self):
+        self.summary_tsv.resolve().parent.mkdir(parents=True, exist_ok=True)
+
+        if self.summary_tsv.exists():
+            with open(self.summary_tsv, "r", encoding="utf-8") as f:
+                last_row = f.readlines()[-1]
+        else:
+            self.write_summary_tsv(
+                [
+                    'method',
+                    '# candidates',
+                    'predecessor change',
+                    'current model criterion',
+                    'current model',
+                    'candidate changes',
+                ]
+            )
+
+    @classmethod
+    def read_arguments_from_yaml_dict(cls, yaml_dict):
+        kwargs = copy.deepcopy(yaml_dict)
+
+        predecessor_model = None
+        if (
+            predecessor_model_yaml := kwargs.pop(PREDECESSOR_MODEL, None)
+        ) is not None:
+
+            predecessor_model = Model.from_yaml(predecessor_model_yaml)
+
+        return {
+            **kwargs,
+            PREDECESSOR_MODEL: predecessor_model,
+        }
+
     def is_plausible(self, model: Model) -> bool:
         """Determine whether a candidate model is plausible.
 
@@ -93,7 +188,7 @@ def accept(
         self,
         model: Model,
         distance: Union[None, float, int],
-        #keep_others: bool = True,
+        # keep_others: bool = True,
     ) -> None:
         """Add a candidate model to the candidate space.
 
@@ -135,6 +230,15 @@ def exclude(
         else:
             self.exclusions.append(model.get_hash())
 
+    def exclude_hashes(self, hashes: List[str]) -> None:
+        """Exclude models from future consideration, by hash.
+
+        Args:
+            hashes:
+                The model hashes that will be excluded
+        """
+        self.exclusions.extend(hashes)
+
     def excluded(
         self,
         model: Model,
@@ -178,13 +282,15 @@ def consider(self, model: Union[Model, None]) -> bool:
         if self.limit.reached():
             return False
         if self.excluded(model):
-            warnings.warn(f'Model has been previously excluded from the candidate space so is skipped here. Model subspace ID: {model.model_subspace_id}. Parameterization: {model.parameters}',  RuntimeWarning)
+            warnings.warn(
+                f'Model has been previously excluded from the candidate space so is skipped here. Model subspace ID: {model.model_subspace_id}. Parameterization: {model.parameters}',
+                RuntimeWarning,
+            )
             return True
         if not self.is_plausible(model):
             return True
         if not self._consider_method(model):
             return True
-
         self.accept(model, distance=self.distance(model))
         return not self.limit.reached()
 
@@ -193,18 +299,22 @@ def reset_accepted(self) -> None:
         self.models = []
         self.distances = []
 
-    def set_predecessor_model(self, predecessor_model: Union[Model, str, None]):
+    def set_predecessor_model(
+        self, predecessor_model: Union[Model, str, None]
+    ):
         self.predecessor_model = predecessor_model
         if (
             self.predecessor_model == VIRTUAL_INITIAL_MODEL
             and self.method not in VIRTUAL_INITIAL_MODEL_METHODS
         ):
-            raise ValueError(f'A virtual initial model was requested for a method ({self.method}) that does not support them.')  # noqa: E501
+            raise ValueError(
+                f'A virtual initial model was requested for a method ({self.method}) that does not support them.'
+            )
 
     def get_predecessor_model(self):
         return self.predecessor_model
 
-    def set_exclusions(self, exclusions: Union[List[Any], None]):
+    def set_exclusions(self, exclusions: Union[List[str], None]):
         # TODO change to List[str] for hashes?
         self.exclusions = exclusions
         if self.exclusions is None:
@@ -220,11 +330,27 @@ def set_limit(self, limit: TYPE_LIMIT = None):
     def get_limit(self):
         return self.limit.get_limit()
 
+    def wrap_search_subspaces(self, search_subspaces: Callable[[], None]):
+        """Decorate the subspace searches of a model space.
+
+        Used by candidate spaces to perform changes that alter the search.
+        See `BidirectionalCandidateSpace` for an example, where it's used to switch directions.
+
+        Args:
+            search_subspaces:
+                The method that searches the subspaces of a model space.
+        """
+
+        def wrapper():
+            search_subspaces()
+
+        return wrapper
+
     def reset(
         self,
         predecessor_model: Optional[Union[Model, str, None]] = None,
         # FIXME change `Any` to some `TYPE_MODEL_HASH` (e.g. union of str/int/float)
-        exclusions: Optional[Union[List[Any], None]] = None,
+        exclusions: Optional[Union[List[str], None]] = None,
         limit: TYPE_LIMIT = None,
     ) -> None:
         """Reset the candidate models, optionally reinitialize with a model.
@@ -271,9 +397,8 @@ def distances_in_estimated_parameters(
             model0 = self.predecessor_model
         model1 = model
 
-        if (
-            model0 != VIRTUAL_INITIAL_MODEL
-            and not model1.petab_yaml.samefile(model0.petab_yaml)
+        if model0 != VIRTUAL_INITIAL_MODEL and not model1.petab_yaml.samefile(
+            model0.petab_yaml
         ):
             raise NotImplementedError(
                 'Computation of distances between different PEtab problems is '
@@ -300,9 +425,13 @@ def distances_in_estimated_parameters(
                 )
         else:
             parameter_ids = list(model0.petab_parameters)
-            parameters0 = np.array(model0.get_parameter_values(parameter_ids=parameter_ids))
+            parameters0 = np.array(
+                model0.get_parameter_values(parameter_ids=parameter_ids)
+            )
 
-        parameters1 = np.array(model1.get_parameter_values(parameter_ids=parameter_ids))
+        parameters1 = np.array(
+            model1.get_parameter_values(parameter_ids=parameter_ids)
+        )
 
         # TODO change to some numpy elementwise operation
         estimated0 = np.array([p == ESTIMATE for p in parameters0]).astype(int)
@@ -320,10 +449,25 @@ def distances_in_estimated_parameters(
         distances = {
             'l1': l1,
             'size': size,
-           }
+        }
 
         return distances
 
+    def update_after_calibration(
+        self,
+        *args,
+        **kwargs,
+    ):
+        """Do work in the candidate space after calibration.
+
+        For example, this is used by the `FamosCandidateSpace` to switch
+        methods.
+
+        Different candidate spaces require different arguments. All arguments
+        are here, to ensure candidate spaces can be switched easily and still
+        receive sufficient arguments.
+        """
+        pass
 
 
 class ForwardCandidateSpace(CandidateSpace):
@@ -332,7 +476,11 @@ class ForwardCandidateSpace(CandidateSpace):
     Attributes:
         direction:
             `1` for the forward method, `-1` for the backward method.
+        max_steps:
+            Maximum number of steps forward in a single iteration of forward selection.
+            Defaults to no maximum (`None`).
     """
+
     method = Method.FORWARD
     direction = 1
 
@@ -340,28 +488,31 @@ def __init__(
         self,
         *args,
         predecessor_model: Optional[Union[Model, str]] = None,
+        max_steps: int = None,
         **kwargs,
     ):
         # Although `VIRTUAL_INITIAL_MODEL` is `str` and can be used as a default
         # argument, `None` may be passed by other packages, so the default value
         # is handled here instead.
+        self.max_steps = max_steps
         if predecessor_model is None:
             predecessor_model = VIRTUAL_INITIAL_MODEL
         super().__init__(*args, predecessor_model=predecessor_model, **kwargs)
 
     def is_plausible(self, model: Model) -> bool:
-        distances = \
-            self.distances_in_estimated_parameters(model)
-        unsigned_size = self.direction * distances['size']
+        distances = self.distances_in_estimated_parameters(model)
+        n_steps = self.direction * distances['size']
+
+        if self.max_steps is not None and n_steps > self.max_steps:
+            raise StopIteration(
+                f"Maximal number of steps for method {self.method} exceeded. Stop sending candidate models."
+            )
+
         # A model is plausible if the number of estimated parameters strictly
         # increases (or decreases, if `self.direction == -1`), and no
         # previously estimated parameters become fixed.
-        if (
-            self.predecessor_model == VIRTUAL_INITIAL_MODEL
-            or (
-                unsigned_size > 0
-                and distances['l1'] == unsigned_size
-            )
+        if self.predecessor_model == VIRTUAL_INITIAL_MODEL or (
+            n_steps > 0 and distances['l1'] == n_steps
         ):
             return True
         return False
@@ -369,8 +520,7 @@ def is_plausible(self, model: Model) -> bool:
     def distance(self, model: Model) -> int:
         # TODO calculated here and `is_plausible`. Rewrite to only calculate
         #      once?
-        distances = \
-            self.distances_in_estimated_parameters(model)
+        distances = self.distances_in_estimated_parameters(model)
         return distances['l1']
 
     def _consider_method(self, model) -> bool:
@@ -382,8 +532,8 @@ def _consider_method(self, model) -> bool:
         if self.distances:
             distance0 = one(set(self.distances))
             # TODO store each or just one?
-            #distance0 = one(self.distances)
-        #breakpoint()
+            # distance0 = one(self.distances)
+        # breakpoint()
 
         # Only keep the best model(s).
         if distance > distance0:
@@ -395,29 +545,779 @@ def _consider_method(self, model) -> bool:
 
 class BackwardCandidateSpace(ForwardCandidateSpace):
     """The backward method class."""
+
     method = Method.BACKWARD
     direction = -1
 
 
-class LateralCandidateSpace(ForwardCandidateSpace):
+class BidirectionalCandidateSpace(ForwardCandidateSpace):
+    """The bidirectional method class.
+
+    Attributes:
+        method_history:
+            The history of models that were found at each search.
+            A list of dictionaries, where each dictionary contains keys for the `METHOD`
+            and the list of `MODELS`.
+    """
+
+    # TODO refactor method to inherit from governing_method if not specified
+    #      by constructor argument -- remove from here.
+    method = Method.BIDIRECTIONAL
+    governing_method = Method.BIDIRECTIONAL
+    retry_model_space_search_if_no_models = True
+
+    def __init__(
+        self,
+        *args,
+        initial_method: Method = Method.FORWARD,
+        **kwargs,
+    ):
+        super().__init__(*args, **kwargs)
+
+        # FIXME cannot access from CLI
+        # FIXME probably fine to replace `self.initial_method`
+        #       with `self.method` here. i.e.:
+        #       1. change `method` to `Method.FORWARD
+        #       2. change signature to `initial_method: Method = None`
+        #       3. change code here to `if initial_method is not None: self.method = initial_method`
+        self.initial_method = initial_method
+
+        self.method_history: List[Dict[str, Union[Method, List[Model]]]] = []
+
+    def update_method(self, method: Method):
+        if method == Method.FORWARD:
+            self.direction = 1
+        elif method == Method.BACKWARD:
+            self.direction = -1
+        else:
+            raise NotImplementedError(
+                f'Bidirectional direction must be either `Method.FORWARD` or `Method.BACKWARD`, not {method}.'
+            )
+
+        self.method = method
+
+    def switch_method(self):
+        if self.method == Method.FORWARD:
+            method = Method.BACKWARD
+        elif self.method == Method.BACKWARD:
+            method = Method.FORWARD
+
+        self.update_method(method=method)
+
+    def setup_before_model_subspaces_search(self):
+        # If previous search found no models, then switch method.
+        previous_search = (
+            None if not self.method_history else self.method_history[-1]
+        )
+        if previous_search is None:
+            self.update_method(self.initial_method)
+            return
+
+        self.update_method(previous_search[METHOD])
+        if not previous_search[MODELS]:
+            self.switch_method()
+            self.retry_model_space_search_if_no_models = False
+
+    def setup_after_model_subspaces_search(self):
+        current_search = {
+            METHOD: self.method,
+            MODELS: self.models,
+        }
+        self.method_history.append(current_search)
+        self.method = self.governing_method
+
+    def wrap_search_subspaces(self, search_subspaces):
+        def wrapper():
+            def iterate():
+                self.setup_before_model_subspaces_search()
+                search_subspaces()
+                self.setup_after_model_subspaces_search()
+
+            # Repeat until models are found or switching doesn't help.
+            iterate()
+            while (
+                not self.models and self.retry_model_space_search_if_no_models
+            ):
+                iterate()
+
+            # Reset flag for next time.
+            self.retry_model_space_search_if_no_models = True
+
+        return wrapper
+
+
+class FamosCandidateSpace(CandidateSpace):
+    """The FAMoS method class.
+
+    This candidate space implements and extends the original FAMoS
+    algorithm (doi: 10.1371/journal.pcbi.1007230).
+
+    Attributes:
+        critical_parameter_sets:
+            A list of lists, where each inner list contains parameter IDs.
+            All models must estimate at least 1 parameter from each critical
+            parameter set.
+        swap_parameter_sets:
+            A list of lists, where each inner list contains parameter IDs.
+            The lateral moves in FAMoS are constrained to be between parameters that
+            exist in the same swap parameter set.
+        method_scheme:
+            A dictionary that specifies how to switch between methods when
+            the current method doesn't produce a better model.
+            Keys are `n`-tuples that described a pattern of length `n`
+            methods. Values are methods. If the previous methods match the
+            pattern in the key, then the method in the value will be used next.
+            The order of the dictionary is important: only the first matched
+            pattern will be used.
+            Defaults to the method scheme described in the original FAMoS
+            publication.
+        n_reattempts:
+            Integer. The total number of times that a jump-to-most-distance action
+            will be performed, triggered whenever the model selection would
+            normally terminate. Defaults to no reattempts (`0`).
+        consecutive_laterals:
+            Boolean. If `True`, the method will continue performing lateral moves
+            while they produce better models. Otherwise, the method scheme will
+            be applied after one lateral move.
+    """
+
+    method = Method.FAMOS
+    default_method_scheme = {
+        (Method.BACKWARD, Method.FORWARD): Method.LATERAL,
+        (Method.FORWARD, Method.BACKWARD): Method.LATERAL,
+        (Method.BACKWARD, Method.LATERAL): None,
+        (Method.FORWARD, Method.LATERAL): None,
+        (Method.FORWARD,): Method.BACKWARD,
+        (Method.BACKWARD,): Method.FORWARD,
+        (Method.LATERAL,): Method.FORWARD,
+        (Method.MOST_DISTANT,): Method.FORWARD,
+        None: Method.FORWARD,
+    }
+
+    def __init__(
+        self,
+        *args,
+        predecessor_model: Optional[Union[Model, str, None]] = None,
+        critical_parameter_sets: List = [],
+        swap_parameter_sets: List = [],
+        method_scheme: Dict[tuple, str] = None,
+        n_reattempts: int = 0,
+        consecutive_laterals: bool = False,
+        **kwargs,
+    ):
+        self.critical_parameter_sets = critical_parameter_sets
+        self.swap_parameter_sets = swap_parameter_sets
+
+        self.method_scheme = method_scheme
+        if method_scheme is None:
+            self.method_scheme = self.default_method_scheme
+
+        # FIXME remove and use `self.method` everywhere in the constructor
+        #       instead -- not required in other methods anymore
+        self.initial_method = self.method_scheme[None]
+        self.method = self.initial_method
+        self.method_history = [self.initial_method]
+
+        if predecessor_model is None:
+            predecessor_model = VIRTUAL_INITIAL_MODEL
+
+        if (
+            predecessor_model == VIRTUAL_INITIAL_MODEL
+            and critical_parameter_sets
+        ) or (
+            predecessor_model != VIRTUAL_INITIAL_MODEL
+            and not self.check_critical(predecessor_model)
+        ):
+            raise ValueError(
+                f'Provided predecessor model {predecessor_model.parameters} does not contain necessary critical parameters {self.critical_parameter_sets}. Provide a valid predecessor model.'
+            )
+
+        if (
+            predecessor_model == VIRTUAL_INITIAL_MODEL
+            and self.initial_method not in VIRTUAL_INITIAL_MODEL_METHODS
+        ):
+            raise ValueError(
+                f"The initial method {self.initial_method} does not support the `VIRTUAL_INITIAL_MODEL` as its predecessor_model."
+            )
+
+        # FIXME remove `None` from the resulting `inner_methods` set?
+        inner_methods = set.union(
+            *[
+                set(
+                    [
+                        *(
+                            method_pattern
+                            if method_pattern is not None
+                            else (None,)
+                        ),
+                        next_method,
+                    ]
+                )
+                for method_pattern, next_method in self.method_scheme.items()
+            ]
+        )
+        if Method.LATERAL in inner_methods and not self.swap_parameter_sets:
+            raise ValueError(
+                f"Use of the lateral method with FAMoS requires `swap_parameter_sets`."
+            )
+
+        for method in inner_methods:
+            if method is not None and method not in [
+                Method.FORWARD,
+                Method.LATERAL,
+                Method.BACKWARD,
+                Method.MOST_DISTANT,
+            ]:
+                raise NotImplementedError(
+                    f'Methods FAMoS can swap to are `Method.FORWARD`, `Method.BACKWARD` and `Method.LATERAL`, not {method}. \
+                    Check if the method_scheme scheme provided is correct.'
+                )
+
+        self.inner_candidate_spaces = {
+            Method.FORWARD: ForwardCandidateSpace(
+                *args,
+                predecessor_model=predecessor_model,
+                **kwargs,
+            ),
+            Method.BACKWARD: BackwardCandidateSpace(
+                *args,
+                predecessor_model=predecessor_model,
+                **kwargs,
+            ),
+            Method.LATERAL: LateralCandidateSpace(
+                *args,
+                predecessor_model=(
+                    predecessor_model
+                    if predecessor_model != VIRTUAL_INITIAL_MODEL
+                    else None
+                ),
+                max_steps=1,
+                **kwargs,
+            ),
+        }
+        self.inner_candidate_space = self.inner_candidate_spaces[
+            self.initial_method
+        ]
+
+        super().__init__(*args, predecessor_model=predecessor_model, **kwargs)
+
+        self.governing_method = Method.FAMOS
+
+        self.n_reattempts = n_reattempts
+
+        self.consecutive_laterals = consecutive_laterals
+        if (
+            not self.consecutive_laterals
+            and (Method.LATERAL,) not in self.method_scheme
+        ):
+            raise ValueError(
+                "Please provide a method to switch to after a lateral search, if not enabling the `consecutive_laterals` option."
+            )
+
+        if self.n_reattempts:
+            # TODO make so max_number can be specified? It cannot in original FAMoS.
+            self.most_distant_max_number = 100
+        else:
+            self.most_distant_max_number = 1
+
+        self.best_models = []
+        self.best_model_of_current_run = predecessor_model
+
+        self.jumped_to_most_distant = False
+        self.swap_done_successfully = False
+
+    @classmethod
+    def read_arguments_from_yaml_dict(cls, yaml_dict) -> dict:
+        kwargs = copy.deepcopy(yaml_dict)
+
+        if (method_scheme_raw := kwargs.pop(METHOD_SCHEME, None)) is not None:
+            method_scheme = {
+                (
+                    tuple(
+                        [
+                            Method(method_str)
+                            for method_str in definition[PREVIOUS_METHODS]
+                        ]
+                    )
+                    if definition[PREVIOUS_METHODS] is not None
+                    else None
+                ): definition[NEXT_METHOD]
+                for definition in method_scheme_raw
+            }
+            kwargs[METHOD_SCHEME] = method_scheme
+
+        return super().read_arguments_from_yaml_dict(yaml_dict=kwargs)
+
+    def update_after_calibration(
+        self,
+        *args,
+        calibrated_models: Dict[str, Model],
+        newly_calibrated_models: Dict[str, Model],
+        criterion: Criterion,
+        **kwargs,
+    ) -> None:
+        """See `CandidateSpace.update_after_calibration`."""
+
+        # super().update_after_calibration(*args, **kwargs)
+
+        # In case we jumped to most distant in the last iteration,
+        # there's no need for an update, so we reset the jumped variable
+        # to False and continue to candidate generation
+        if self.jumped_to_most_distant:
+            self.jumped_to_most_distant = False
+            jumped_to_model = one(newly_calibrated_models.values())
+            self.set_predecessor_model(jumped_to_model)
+            self.best_model_of_current_run = jumped_to_model
+            return False
+
+        if self.update_from_newly_calibrated_models(
+            newly_calibrated_models=newly_calibrated_models,
+            criterion=criterion,
+        ):
+            logging.info("Switching method")
+            self.switch_method(calibrated_models=calibrated_models)
+            self.switch_inner_candidate_space(
+                calibrated_models=calibrated_models,
+            )
+            logging.info(
+                "Method switched to ", self.inner_candidate_space.method
+            )
+
+        self.method_history.append(self.method)
+
+    def update_from_newly_calibrated_models(
+        self,
+        newly_calibrated_models: Dict[str, Model],
+        criterion: Criterion,
+    ) -> bool:
+        """Update the self.best_models with the latest
+        `newly_calibrated_models`
+        and determine if there was a new best model. If so, return
+        True. False otherwise."""
+
+        go_into_switch_method = True
+        for newly_calibrated_model in newly_calibrated_models.values():
+            if (
+                self.best_model_of_current_run == VIRTUAL_INITIAL_MODEL
+                or default_compare(
+                    model0=self.best_model_of_current_run,
+                    model1=newly_calibrated_model,
+                    criterion=criterion,
+                )
+            ):
+                go_into_switch_method = False
+                self.best_model_of_current_run = newly_calibrated_model
+
+            if len(
+                self.best_models
+            ) < self.most_distant_max_number or default_compare(
+                model0=self.best_models[self.most_distant_max_number - 1],
+                model1=newly_calibrated_model,
+                criterion=criterion,
+            ):
+                self.insert_model_into_best_models(
+                    model_to_insert=newly_calibrated_model,
+                    criterion=criterion,
+                )
+
+        self.best_models = self.best_models[: self.most_distant_max_number]
+
+        # When we switch to LATERAL method, we will do only one iteration with this
+        # method. So if we do it succesfully (i.e. that we found a new best model), we
+        # want to switch method. This is why we put go_into_switch_method to True, so
+        # we go into the method switching pipeline
+        if self.method == Method.LATERAL and not self.consecutive_laterals:
+            self.swap_done_successfully = True
+            go_into_switch_method = True
+        return go_into_switch_method
+
+    def insert_model_into_best_models(
+        self, model_to_insert: Model, criterion: Criterion
+    ) -> None:
+        """Inserts a model into the list of best_models which are sorted
+        w.r.t. the criterion specified."""
+        insert_index = bisect.bisect_left(
+            [model.get_criterion(criterion) for model in self.best_models],
+            model_to_insert.get_criterion(criterion),
+        )
+        self.best_models.insert(insert_index, model_to_insert)
+
+    def consider(self, model: Union[Model, None]) -> bool:
+        """Re-define consider of FAMoS to be the consider method
+        of the inner_candidate_space. Update all the attributes
+        changed in the condsider method."""
+
+        if self.limit.reached():
+            return False
+
+        # Check if model contains necessary critical parameters and, if
+        # the current method is swap, check the swap move is contained
+        # in a swap set.
+
+        if self.check_critical(model) and (
+            not self.method == Method.LATERAL or self.check_swap(model)
+        ):
+            return_value = self.inner_candidate_space.consider(model)
+
+            # update the attributes
+            self.models = self.inner_candidate_space.models
+            self.distances = self.inner_candidate_space.distances
+            self.exclusions = self.inner_candidate_space.exclusions
+
+            return return_value
+
+        return True
+
+    def _consider_method(self, model) -> bool:
+        """See `CandidateSpace._consider_method`."""
+        return self.inner_candidate_space._consider_method(model)
+
+    def reset_accepted(self) -> None:
+        """Changing the reset_accepted to reset the
+        inner_candidate_space as well."""
+        super().reset_accepted()
+        self.inner_candidate_space.reset_accepted()
+
+    def set_predecessor_model(
+        self, predecessor_model: Union[Model, str, None]
+    ):
+        """Setting the predecessor model for the
+        inner_candidate_space as well."""
+        super().set_predecessor_model(predecessor_model=predecessor_model)
+        self.inner_candidate_space.set_predecessor_model(
+            predecessor_model=predecessor_model
+        )
+
+    def set_exclusions(self, exclusions: Union[List[str], None]):
+        """Setting the exclusions for the
+        inner_candidate_space as well."""
+        self.exclusions = exclusions
+        self.inner_candidate_space.exclusions = exclusions
+        if self.exclusions is None:
+            self.exclusions = []
+            self.inner_candidate_space.exclusions = []
+
+    def set_limit(self, limit: TYPE_LIMIT = None):
+        """Setting the limit for the
+        inner_candidate_space as well."""
+        if limit is not None:
+            self.limit.set_limit(limit)
+            self.inner_candidate_space.limit.set_limit(limit)
+
+    def is_plausible(self, model: Model) -> bool:
+        if self.check_critical(model):
+            return self.inner_candidate_space.is_plausible(model)
+        return False
+
+    def check_swap(self, model: Model) -> bool:
+        """Check if parameters that are swapped are contained in the
+        same swap parameter set."""
+        if self.method != Method.LATERAL:
+            return True
+
+        predecessor_estimated_parameters_ids = set(
+            self.predecessor_model.get_estimated_parameter_ids_all()
+        )
+        estimated_parameters_ids = set(model.get_estimated_parameter_ids_all())
+
+        swapped_parameters_ids = estimated_parameters_ids.symmetric_difference(
+            predecessor_estimated_parameters_ids
+        )
+
+        for swap_set in self.swap_parameter_sets:
+            if swapped_parameters_ids.issubset(set(swap_set)):
+                return True
+        return False
+
+    def check_critical(self, model: Model) -> bool:
+        """Check if the model contains all necessary critical parameters"""
+
+        estimated_parameters_ids = set(model.get_estimated_parameter_ids_all())
+        for critical_set in self.critical_parameter_sets:
+            if not estimated_parameters_ids.intersection(set(critical_set)):
+                return False
+        return True
+
+    def switch_method(
+        self,
+        calibrated_models: Dict[str, Model],
+    ) -> None:
+        """Switch to the next method with respect to the history
+        of methods used and the switching scheme in self.method_scheme"""
+
+        previous_method = self.method
+        next_method = previous_method
+        logging.info("SWITCHING", self.method_history)
+        # If last method was LATERAL and we have made a succesfull swap
+        # (a better model was found) then go back to FORWARD. Else do the
+        # usual method swapping scheme.
+        if self.swap_done_successfully:
+            next_method = self.method_scheme[(Method.LATERAL,)]
+        else:
+            # iterate through the method_scheme dictionary to see which method to switch to
+            for previous_methods in self.method_scheme:
+                if previous_methods is not None and previous_methods == tuple(
+                    self.method_history[-len(previous_methods) :]
+                ):
+                    next_method = self.method_scheme[previous_methods]
+                    # if found a switch just break (choosing first good switch)
+                    break
+
+        # raise error if the method didn't change
+        if next_method == previous_method:
+            raise ValueError(
+                "Method didn't switch when it had to. "
+                "The `method_scheme` provided is not sufficient. "
+                "Please provide a correct method_scheme scheme. "
+                f"Method history: `{self.method_history}`. "
+            )
+
+        # Terminate if next method is `None`
+        if next_method is None:
+            if self.n_reattempts:
+                self.jump_to_most_distant(calibrated_models=calibrated_models)
+                return
+            raise StopIteration(
+                f"The next method is {next_method}. The search is terminating."
+            )
+
+        # If we try to switch to SWAP method but it's not available (no crit or swap parameter groups)
+        if (
+            next_method == Method.LATERAL
+            and not [
+                critical_set
+                for critical_set in self.critical_parameter_sets
+                if len(critical_set) > 1
+            ]
+            and not self.swap_parameter_sets
+        ):
+            if self.n_reattempts:
+                self.jump_to_most_distant(calibrated_models=calibrated_models)
+                return
+            raise StopIteration(
+                f"The next method is {next_method}, but there are no critical or swap parameters sets. Terminating."
+            )
+        if previous_method == Method.LATERAL:
+            self.swap_done_successfully = False
+        self.update_method(method=next_method)
+
+    def update_method(self, method: Method):
+        """Update self.method to the method."""
+
+        self.method = method
+
+    def switch_inner_candidate_space(
+        self,
+        calibrated_models: Dict[str, Model],
+    ):
+        """Switch self.inner_candidate_space to the candidate space of
+        the current self.method."""
+
+        # if self.method != Method.MOST_DISTANT:
+        self.inner_candidate_space = self.inner_candidate_spaces[self.method]
+        # reset the next inner candidate space with the current history of all
+        # calibrated models
+        self.inner_candidate_space.reset(
+            predecessor_model=self.predecessor_model,
+            exclusions=calibrated_models,
+        )
+
+    def jump_to_most_distant(
+        self,
+        calibrated_models: Dict[str, Model],
+    ):
+        """Jump to most distant model with respect to the history of all
+        calibrated models."""
+
+        predecessor_model = self.get_most_distant(
+            calibrated_models=calibrated_models,
+        )
+
+        logging.info("JUMPING: ", predecessor_model.parameters)
+
+        # if model not appropriate make it so by adding the first
+        # critical parameter from each critical parameter set
+        if not self.check_critical(predecessor_model):
+            for critical_set in self.critical_parameter_sets:
+                predecessor_model.parameters[critical_set[0]] = ESTIMATE
+
+        # self.update_method(self.initial_method)
+        self.update_method(Method.MOST_DISTANT)
+
+        self.n_reattempts -= 1
+        self.jumped_to_most_distant = True
+
+        # FIXME rename here `predecessor_model` to `most_distant_model`
+        # self.predecessor_model = None
+        self.set_predecessor_model(None)
+        self.best_model_of_current_run = None
+        self.models = [predecessor_model]
+
+        self.write_summary_tsv("Jumped to the most distant model.")
+        self.update_method(self.method_scheme[(Method.MOST_DISTANT,)])
+
+    # TODO Fix for non-famos model subspaces. FAMOS easy because of only 0;ESTIMATE
+    def get_most_distant(
+        self,
+        calibrated_models: Dict[str, Model],
+    ) -> Model:
+        """
+        Get most distant model to all the checked models. We take models from the
+        sorted list of best models (self.best_models) and construct complements of
+        these models. For all these complements we compute the distance in number of
+        different estimated parameters to all models from history. For each complement
+        we take the minimum of these distances as it's distance to history. Then we
+        choose the complement model with the maximal distance to history.
+        TODO:
+        Next we check if this model is contained in any subspace. If so we choose it.
+        If not we choose the model in a subspace that has least distance to this
+        complement model.
+        """
+
+        most_distance = 0
+        most_distant_indices = []
+
+        parameter_ids = self.best_models[0].petab_parameters
+
+        for model in self.best_models:
+            model_estimated_parameters = np.array(
+                [
+                    p == ESTIMATE
+                    for p in model.get_parameter_values(
+                        parameter_ids=parameter_ids
+                    )
+                ]
+            ).astype(int)
+            complement_parameters = 1 - model_estimated_parameters
+            # initialize the least distance to the maximal possible value of it
+            complement_least_distance = len(complement_parameters)
+            # get the complement least distance
+            for calibrated_model in calibrated_models.values():
+                calibrated_model_estimated_parameters = np.array(
+                    [
+                        p == ESTIMATE
+                        for p in calibrated_model.get_parameter_values(
+                            parameter_ids=parameter_ids
+                        )
+                    ]
+                ).astype(int)
+
+                difference = (
+                    calibrated_model_estimated_parameters
+                    - complement_parameters
+                )
+                l1_distance = np.abs(difference).sum()
+
+                if l1_distance < complement_least_distance:
+                    complement_least_distance = l1_distance
+
+            # if the complement is further away than the current best one
+            # then save it as the new most distant away one
+            if complement_least_distance > most_distance:
+                most_distance = complement_least_distance
+                most_distant_indices = complement_parameters
+        if len(most_distant_indices) == 0:
+            raise StopIteration("No most_distant model found. Terminating")
+
+        most_distant_parameter_values = [
+            str(index).replace('1', ESTIMATE) for index in most_distant_indices
+        ]
+
+        most_distant_parameters = {
+            parameter_id: index
+            for parameter_id, index in zip(
+                parameter_ids, most_distant_parameter_values
+            )
+        }
+
+        most_distant_model = Model(
+            petab_yaml=model.petab_yaml,
+            model_subspace_id=model.model_subspace_id,
+            model_subspace_indices=most_distant_indices,
+            parameters=most_distant_parameters,
+        )
+
+        return most_distant_model
+
+    def wrap_search_subspaces(self, search_subspaces):
+        def wrapper():
+            search_subspaces(only_one_subspace=True)
+
+        return wrapper
+
+
+# TODO rewrite so BidirectionalCandidateSpace inherits from ForwardAndBackwardCandidateSpace
+#      instead
+class ForwardAndBackwardCandidateSpace(BidirectionalCandidateSpace):
+    method = Method.FORWARD_AND_BACKWARD
+    governing_method = Method.FORWARD_AND_BACKWARD
+    retry_model_space_search_if_no_models = False
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs, initial_method=None)
+
+    def wrap_search_subspaces(self, search_subspaces):
+        def wrapper():
+            for method in [Method.FORWARD, Method.BACKWARD]:
+                self.update_method(method=method)
+                search_subspaces()
+                self.setup_after_model_subspaces_search()
+
+        return wrapper
+
+    # Disable unused interface
+    setup_before_model_subspaces_search = None
+    switch_method = None
+
+    def setup_after_model_space_search(self):
+        pass
+
+
+class LateralCandidateSpace(CandidateSpace):
     """Find models with the same number of estimated parameters."""
+
     method = Method.LATERAL
 
     def __init__(
         self,
         *args,
-        predecessor_model: Union[Model, str],
+        predecessor_model: Union[Model, None],
+        max_steps: int = None,
         **kwargs,
     ):
-        super().__init__(*args, predecessor_model=predecessor_model, **kwargs)
+        """
+        Additional args:
+            max_number_of_steps:
+                Maximal allowed number of swap moves. If 0 then there is no maximum.
+        """
+        super().__init__(
+            *args,
+            predecessor_model=predecessor_model,
+            **kwargs,
+        )
+        self.max_steps = max_steps
 
     def is_plausible(self, model: Model) -> bool:
+        if self.predecessor_model is None:
+            raise ValueError(
+                f"The predecessor_model is still None. Provide an appropriate predecessor_model"
+            )
+
         distances = self.distances_in_estimated_parameters(model)
+
+        # If max_number_of_steps is non-zero and the number of steps made is
+        # larger then move is not plausible.
+        if self.max_steps is not None and distances['l1'] > 2 * self.max_steps:
+            raise StopIteration(
+                f"Maximal number of steps for method {self.method} exceeded. Stop sending candidate models."
+            )
+
         # A model is plausible if the number of estimated parameters remains
         # the same, but some estimated parameters have become fixed and vice
         # versa.
         if (
-            distances['size'] == 0 and
+            distances['size'] == 0
+            and
             # distances['size'] == 0 implies L1 % 2 == 0.
             # FIXME here and elsewhere, deal with models that are equal
             #       except for the values of their fixed parameters.
@@ -432,10 +1332,11 @@ def _consider_method(self, model):
 
 class BruteForceCandidateSpace(CandidateSpace):
     """The brute-force method class."""
+
     method = Method.BRUTE_FORCE
 
     def __init__(self, *args, **kwargs):
-        #if args or kwargs:
+        # if args or kwargs:
         #    # FIXME remove?
         #    # FIXME at least support limit
         #    warnings.warn(
@@ -451,8 +1352,11 @@ def _consider_method(self, model):
 candidate_space_classes = [
     ForwardCandidateSpace,
     BackwardCandidateSpace,
+    BidirectionalCandidateSpace,
     LateralCandidateSpace,
     BruteForceCandidateSpace,
+    ForwardAndBackwardCandidateSpace,
+    FamosCandidateSpace,
 ]
 
 
@@ -470,7 +1374,9 @@ def method_to_candidate_space_class(method: Method) -> str:
     for candidate_space_class in candidate_space_classes:
         if candidate_space_class.method == method:
             return candidate_space_class
-    raise NotImplementedError(f'The provided method name {method} does not correspond to an implemented candidate space.')  # noqa: E501
+    raise NotImplementedError(
+        f'The provided method name {method} does not correspond to an implemented candidate space.'
+    )
 
 
 '''
diff --git a/petab_select/cli.py b/petab_select/cli.py
index ce8b5597..96b281a5 100644
--- a/petab_select/cli.py
+++ b/petab_select/cli.py
@@ -1,35 +1,65 @@
 """The PEtab Select command-line interface."""
-from pathlib import Path
-from typing import List
 import warnings
-import yaml
+from pathlib import Path
+from typing import Any, Dict, List
 
 import click
 import dill
 import numpy as np
 import pandas as pd
+import yaml
+from more_itertools import one
 
 from . import ui
-from .constants import (
-    INITIAL_MODEL_METHODS,
-    PETAB_YAML,
-)
-from .candidate_space import (
-    BackwardCandidateSpace,
-    BruteForceCandidateSpace,
-    CandidateSpace,
-    ForwardCandidateSpace,
-    LateralCandidateSpace,
-    method_to_candidate_space_class,
-)
-from .model import (
-    Model,
-    models_from_yaml_list,
-    models_to_yaml_list,
-)
+from .candidate_space import CandidateSpace, method_to_candidate_space_class
+from .constants import INITIAL_MODEL_METHODS, PETAB_YAML
+from .model import Model, models_from_yaml_list, models_to_yaml_list
 from .problem import Problem
 
 
+def read_state(filename: str) -> Dict[str, Any]:
+    with open(filename, 'rb') as f:
+        state = dill.load(f)
+
+    state['problem'] = dill.loads(state['problem'])
+    state['candidate_space'] = dill.loads(state['candidate_space'])
+
+    return state
+
+
+def write_state(
+    state: Dict[str, Any],
+    filename: str,
+) -> Dict[str, Any]:
+    with open(filename, 'wb') as f:
+        dill.dump(state, f)
+
+
+def get_state(
+    problem: Problem,
+    candidate_space: CandidateSpace,
+) -> Dict[str, Any]:
+    state = {
+        'problem': dill.dumps(problem),
+        'candidate_space': dill.dumps(candidate_space),
+    }
+    return state
+
+
+def check_state_compatibility(
+    state: Dict[str, Any],
+    problem: Problem,
+    candidate_space: CandidateSpace,
+):
+    # if state['problem'].method != problem.method:
+    #    warnings.warn(
+    #        'The method in the problem loaded from the state does not match '
+    #        'the method specified in either the PEtab Select YAML file or '
+    #        'specified explicitly with this command.'
+    #    )
+    pass
+
+
 @click.group()
 def cli():
     pass
@@ -37,15 +67,15 @@ def cli():
 
 @cli.command("candidates")
 @click.option(
-    '--yaml',
-    '-y',
-    'yaml_',
+    '--problem',
+    '-p',
+    'problem_yaml',
     help='The PEtab Select YAML problem file.',
 )
 @click.option(
     '--state',
     '-s',
-    'state',
+    'state_filename',
     type=str,
     help='The file that stores the state.',
 )
@@ -62,33 +92,34 @@ def cli():
     'method',
     type=str,
     default=None,
-    help='The method used to identify the candidate models. Defaults to the method in the problem YAML.',  # noqa: E501
+    help='The method used to identify the candidate models. Defaults to the method in the problem YAML.',
 )
 @click.option(
-    '--initial',
-    '-i',
-    'initial',
+    '--previous-predecessor-model',
+    '-P',
+    'previous_predecessor_model_yaml',
     type=str,
     default=None,
-    help='(Optional) The initial model used in the candidate model search.',
+    help='(Optional) The predecessor model used in the previous iteration of model selection.',
 )
 @click.option(
-    '--predecessor',
-    '-p',
-    'predecessor',
+    '--calibrated-models',
+    '-C',
+    'calibrated_models_yamls',
     type=str,
+    multiple=True,
     default=None,
-    help='(Optional) The predecessor model used in the candidate model search.',
+    help='(Optional) Models that have been calibrated.',
 )
 @click.option(
-    '--best',
-    '-b',
-    'best',
+    '--newly-calibrated-models',
+    '-N',
+    'newly_calibrated_models_yamls',
     type=str,
+    multiple=True,
     default=None,
     help=(
-        '(Optional) Use the best model as the predecessor model, from a collection of '
-        'calibrated models.'
+        '(Optional) Models that were calibrated in the most recent iteration.'
     ),
 )
 @click.option(
@@ -118,7 +149,7 @@ def cli():
     help='Whether to output paths relative to the output file.',
 )
 @click.option(
-    '--excluded-model-file',
+    '--excluded-models',
     '-e',
     'excluded_model_files',
     type=str,
@@ -127,7 +158,7 @@ def cli():
     help='Exclude models in this file.',
 )
 @click.option(
-    '--excluded-model-hash-file',
+    '--excluded-model-hashes',
     '-E',
     'excluded_model_hash_files',
     type=str,
@@ -136,13 +167,14 @@ def cli():
     help='Exclude model hashes in this file (one model hash per line).',
 )
 def candidates(
-    yaml_: str,
-    state: str,
+    problem_yaml: str,
+    state_filename: str,
     output: str,
     method: str = None,
-    initial: str = None,
-    predecessor: str = None,
-    best: str = None,
+    previous_predecessor_model_yaml: str = None,
+    # best: str = None,
+    calibrated_models_yamls: List[str] = None,
+    newly_calibrated_models_yamls: List[str] = None,
     limit: float = np.inf,
     limit_sent: float = np.inf,
     relative_paths: bool = False,
@@ -154,25 +186,7 @@ def candidates(
     Documentation for arguments can be viewed with
     `petab_select candidates --help`.
     """
-    if initial is not None:
-        warnings.warn(
-            '`initial` (`-i`) is deprecated in favor of `predecessor` (`-p`), '
-            'for internal consistency.',
-            DeprecationWarning,
-        )
-        if predecessor is not None:
-            raise ValueError(
-                'Only one of `initial` (`-i`) or `predecessor` (`-p`) should be set.'
-            )
-        predecessor = initial
-
-    if predecessor is not None and best is not None:
-        raise KeyError(
-            'The `predecessor` (`-p`) and `best` (`-b`) arguments cannot be used '
-            'together, as they both set the predecessor model.'
-        )
-
-    problem = Problem.from_yaml(yaml_)
+    problem = Problem.from_yaml(problem_yaml)
     if method is None:
         method = problem.method
 
@@ -180,13 +194,20 @@ def candidates(
     # the candidate space.
     problem.method = method
 
-    if Path(state).exists():
-        # Load state
-        with open(state, 'rb') as f:
-            problem.set_state(dill.load(f))
+    candidate_space = problem.new_candidate_space(limit=limit)
+
+    # Setup state
+    if not Path(state_filename).exists():
+        Path(state_filename).parent.mkdir(parents=True, exist_ok=True)
     else:
-        # Create the output path for the state
-        Path(state).parent.mkdir(parents=True, exist_ok=True)
+        state = read_state(state_filename)
+        check_state_compatibility(
+            state=state,
+            problem=problem,
+            candidate_space=candidate_space,
+        )
+        problem = state['problem']
+        candidate_space = state['candidate_space']
 
     excluded_models = []
     # TODO seems like default is `()`, not `None`...
@@ -201,17 +222,44 @@ def candidates(
             with open(excluded_model_hash_file, 'r') as f:
                 excluded_model_hashes += f.read().split('\n')
 
-    predecessor_model = None
-    if best is not None:
-        calibrated_models = models_from_yaml_list(best)
-        predecessor_model = problem.get_best(calibrated_models)
+    previous_predecessor_model = candidate_space.predecessor_model
+    if previous_predecessor_model_yaml is not None:
+        previous_predecessor_model = Model.from_yaml(
+            previous_predecessor_model_yaml
+        )
 
-    if predecessor is not None:
-        predecessor_model = Model.from_yaml(predecessor)
+    # FIXME write single methods to take all models from lists of lists of
+    #       models recursively
+    calibrated_models = None
+    if calibrated_models_yamls is not None:
+        calibrated_models = {}
+        for calibrated_models_yaml in calibrated_models_yamls:
+            calibrated_models.update(
+                {
+                    model.get_hash(): model
+                    for model in models_from_yaml_list(calibrated_models_yaml)
+                }
+            )
 
-    candidate_space = ui.candidates(
+    newly_calibrated_models = None
+    if newly_calibrated_models_yamls is not None:
+        newly_calibrated_models = {}
+        for newly_calibrated_models_yaml in newly_calibrated_models_yamls:
+            newly_calibrated_models.update(
+                {
+                    model.get_hash(): model
+                    for model in models_from_yaml_list(
+                        newly_calibrated_models_yaml
+                    )
+                }
+            )
+
+    ui.candidates(
         problem=problem,
-        predecessor_model=predecessor_model,
+        candidate_space=candidate_space,
+        previous_predecessor_model=previous_predecessor_model,
+        calibrated_models=calibrated_models,
+        newly_calibrated_models=newly_calibrated_models,
         limit=limit,
         limit_sent=limit_sent,
         excluded_models=excluded_models,
@@ -219,8 +267,13 @@ def candidates(
     )
 
     # Save state
-    with open(state, 'wb') as f:
-        dill.dump(problem.get_state(), f)
+    write_state(
+        state=get_state(
+            problem=problem,
+            candidate_space=candidate_space,
+        ),
+        filename=state_filename,
+    )
 
     # Save candidates
     models_to_yaml_list(
@@ -232,9 +285,10 @@ def candidates(
 
 @cli.command("model_to_petab")
 @click.option(
-    '--yaml',
-    '-y',
-    'yaml_',
+    '--model',
+    '-m',
+    'models_yamls',
+    multiple=True,
     help='The PEtab Select model YAML file.',
 )
 @click.option(
@@ -245,8 +299,8 @@ def candidates(
     help='The directory where the PEtab files will be output.',
 )
 @click.option(
-    '--model',
-    '-m',
+    '--model_id',
+    '-i',
     'model_id',
     type=str,
     default=None,
@@ -256,7 +310,7 @@ def candidates(
     ),
 )
 def model_to_petab(
-    yaml_: str,
+    models_yamls: List[str],
     output_path: str,
     model_id: str = None,
 ) -> None:
@@ -267,50 +321,41 @@ def model_to_petab(
     Documentation for arguments can be viewed with
     `petab_select model_to_petab --help`.
     """
+    models = []
+    for models_yaml in models_yamls:
+        models.extend(models_from_yaml_list(models_yaml))
+
+    model0 = None
     try:
-        model = Model.from_yaml(yaml_)
-    except ValueError as e1:
-        # There may be multiple models in a single file.
-        if model_id is None:
-            raise ValueError(
-                'There may be multiple models defined in the provided PEtab '
-                'Select model YAML file. If so, please specify the ID of one '
-                'model that you would like to convert to PEtab, or use '
-                '`models2petab` to convert all of them.'
-            )
-        try:
-            models = models_from_yaml_list(yaml_)
-            try:
-                model = one([
-                    model
-                    for model in models
-                    if model.model_id == model_id
-                ])
-            except ValueError:
-                raise ValueError(
-                    'There must be exactly one model with the provided model '
-                    f'ID "{model_id}", in the provided PEtab Select model '
-                    'YAML file. This may not be the case.'
-                )
-        except Exception as e2:
-            # Unknown error.
-            raise Exception([e1, e2])
-
-    if model_id is not None and model.model_id != model_id:
+        model0 = one(models)
+    except:
+        for model in models:
+            if model.model_id == model_id:
+                if model0 is not None:
+                    raise ValueError('There are multiple models with this ID.')
+                model0 = model
+                # TODO could `break` here and remove the above `ValueError`
+                #      and the `model0` logic
+    if model0 is None:
+        raise ValueError('Could not find a model with the specified model ID.')
+
+    if model_id is not None and model0.model_id != model_id:
         raise ValueError(
             'The ID of the model from the YAML file does not match the '
             'specified ID.'
         )
 
-    result = ui.model_to_petab(model, output_path)
+    result = ui.model_to_petab(model0, output_path)
     print(result[PETAB_YAML])
 
 
 @cli.command("models_to_petab")
 @click.option(
-    '--yaml',
-    '-y',
-    'yaml_',
+    '--models',
+    '-m',
+    'models_yamls',
+    type=str,
+    multiple=True,
     help='The PEtab Select model YAML file, containing a list of models.',
 )
 @click.option(
@@ -321,7 +366,7 @@ def model_to_petab(
     help='The directory where the PEtab files will be output. The PEtab files will be stored in a model-specific subdirectory.',
 )
 def models_to_petab(
-    yaml_: str,
+    models_yamls: List[str],
     output_path_prefix: str,
 ) -> None:
     """Create a PEtab problem for each model in a PEtab Select model YAML file.
@@ -335,7 +380,10 @@ def models_to_petab(
     Documentation for arguments can be viewed with
     `petab_select models_to_petab --help`.
     """
-    models = models_from_yaml_list(yaml_)
+    models = []
+    for models_yaml in models_yamls:
+        models.extend(models_from_yaml_list(models_yaml))
+
     model_ids = pd.Series([model.model_id for model in models])
     duplicates = '\n'.join(set(model_ids[model_ids.duplicated()]))
     if duplicates:
@@ -349,25 +397,29 @@ def models_to_petab(
         models,
         output_path_prefix=output_path_prefix,
     )
-    result_string = '\n'.join([
-        '\t'.join([model.model_id, result[PETAB_YAML]])
-        for model, result in zip(models, results)
-    ])
+    result_string = '\n'.join(
+        [
+            '\t'.join([model.model_id, result[PETAB_YAML]])
+            for model, result in zip(models, results)
+        ]
+    )
     print(result_string)
 
 
 @cli.command("best")
 @click.option(
-    '--yaml',
-    '-y',
-    'yaml_',
+    '--problem',
+    '-p',
+    'problem_yaml',
+    type=str,
     help='The PEtab Select YAML problem file.',
 )
 @click.option(
-    '--models_yaml',
+    '--models',
     '-m',
-    'models_yaml',
+    'models_yamls',
     type=str,
+    multiple=True,
     help='A list of calibrated models.',
 )
 @click.option(
@@ -380,7 +432,7 @@ def models_to_petab(
 @click.option(
     '--state',
     '-s',
-    'state',
+    'state_filename',
     type=str,
     default=None,
     help='The file that stores the state.',
@@ -401,10 +453,10 @@ def models_to_petab(
     help='Whether to output paths relative to the output file.',
 )
 def best(
-    yaml_: str,
-    models_yaml: str,
+    problem_yaml: str,
+    models_yamls: List[str],
     output: str,
-    state: str = None,
+    state_filename: str = None,
     criterion: str = None,
     relative_paths: bool = False,
 ) -> None:
@@ -417,14 +469,30 @@ def best(
     if relative_paths:
         paths_relative_to = Path(output).parent
 
-    problem = Problem.from_yaml(yaml_)
-    if state is not None and Path(state).exists():
-        # Load state
-        with open(state, 'rb') as f:
-            problem.set_state(dill.load(f))
+    problem = Problem.from_yaml(problem_yaml)
+
+    if state_filename is not None:
+        state = read_state(state_filename)
+        check_state_compatibility(
+            state=state,
+            problem=problem,
+            candidate_space=None,
+        )
+        set_state(
+            state=state,
+            problem=problem,
+            candidate_space=None,
+        )
 
-    calibrated_models = models_from_yaml_list(models_yaml)
-    best_model = ui.best(problem=problem, models=calibrated_models, criterion=criterion)
+    models = []
+    for models_yaml in models_yamls:
+        models.extend(models_from_yaml_list(models_yaml))
+
+    best_model = ui.best(
+        problem=problem,
+        models=models,
+        criterion=criterion,
+    )
     best_model.to_yaml(output, paths_relative_to=paths_relative_to)
 
 
diff --git a/petab_select/constants.py b/petab_select/constants.py
index 0e91dfd0..25749482 100644
--- a/petab_select/constants.py
+++ b/petab_select/constants.py
@@ -3,7 +3,6 @@
 from pathlib import Path
 from typing import Dict, List, Union
 
-
 # Zero-indexed column/row indices
 MODEL_ID_COLUMN = 0
 PETAB_YAML_COLUMN = 1
@@ -18,19 +17,20 @@
 PETAB_ESTIMATE_FALSE = 0
 PETAB_ESTIMATE_TRUE = 1
 
-#TYPING_PATH = Union[str, Path]
+# TYPING_PATH = Union[str, Path]
 TYPE_PATH = Union[str, Path]
 
 # Model space file columns
 # TODO ensure none of these occur twice in the column header (this would
 #      suggest that a parameter has a conflicting name)
-#MODEL_ID = 'modelId'  # TODO already defined, reorganize constants
-#YAML = 'YAML'  # FIXME
+# MODEL_ID = 'modelId'  # TODO already defined, reorganize constants
+# YAML = 'YAML'  # FIXME
 MODEL_ID = 'model_id'
 MODEL_SUBSPACE_ID = 'model_subspace_id'
 MODEL_SUBSPACE_INDICES = 'model_subspace_indices'
 MODEL_CODE = 'model_code'
 MODEL_HASH = 'model_hash'
+MODEL_HASHES = 'model_hashes'
 # If `predecessor_model_hash` is defined for a model, it is the ID of the model that the
 # current model was/is to be compared to. This is part of the result and is
 # only (optionally) set by the PEtab calibration tool. It is not defined by the
@@ -42,19 +42,19 @@
 SBML = 'sbml'
 HASH = 'hash'
 
-#MODEL_SPACE_FILE_NON_PARAMETER_COLUMNS = [MODEL_ID, PETAB_YAML]
+# MODEL_SPACE_FILE_NON_PARAMETER_COLUMNS = [MODEL_ID, PETAB_YAML]
 MODEL_SPACE_FILE_NON_PARAMETER_COLUMNS = [MODEL_SUBSPACE_ID, PETAB_YAML]
 
-#COMPARED_MODEL_ID = 'compared_'+MODEL_ID
+# COMPARED_MODEL_ID = 'compared_'+MODEL_ID
 YAML_FILENAME = 'yaml'
 
-#FORWARD = 'forward'
-#BACKWARD = 'backward'
-#BIDIRECTIONAL = 'bidirectional'
-#LATERAL = 'lateral'
+# FORWARD = 'forward'
+# BACKWARD = 'backward'
+# BIDIRECTIONAL = 'bidirectional'
+# LATERAL = 'lateral'
 
 
-#DISTANCES = {
+# DISTANCES = {
 #    FORWARD: {
 #        'l1': 1,
 #        'size': 1,
@@ -67,26 +67,26 @@
 #        'l1': 2,
 #        'size': 0,
 #    },
-#}
+# }
 
 CRITERIA = 'criteria'
 # FIXME remove, change all uses to Enum below
-#AIC = 'AIC'
-#AICC = 'AICc'
-#BIC = 'BIC'
-#AKAIKE_INFORMATION_CRITERION = AIC
-#CORRECTED_AKAIKE_INFORMATION_CRITERION = AICC
-#BAYESIAN_INFORMATION_CRITERION = BIC
-#LH = 'LH'
-#LLH = 'LLH'
-#NLLH = 'NLLH'
-#LIKELIHOOD = LH
-#LOG_LIKELIHOOD = LLH
-#NEGATIVE_LOG_LIKELIHOOD = NLLH
+# AIC = 'AIC'
+# AICC = 'AICc'
+# BIC = 'BIC'
+# AKAIKE_INFORMATION_CRITERION = AIC
+# CORRECTED_AKAIKE_INFORMATION_CRITERION = AICC
+# BAYESIAN_INFORMATION_CRITERION = BIC
+# LH = 'LH'
+# LLH = 'LLH'
+# NLLH = 'NLLH'
+# LIKELIHOOD = LH
+# LOG_LIKELIHOOD = LLH
+# NEGATIVE_LOG_LIKELIHOOD = NLLH
 
 
 PARAMETERS = 'parameters'
-#PARAMETER_ESTIMATE = 'parameter_estimate'
+# PARAMETER_ESTIMATE = 'parameter_estimate'
 ESTIMATED_PARAMETERS = 'estimated_parameters'
 
 CRITERION = 'criterion'
@@ -94,7 +94,14 @@
 VERSION = 'version'
 MODEL_SPACE_FILES = 'model_space_files'
 
+CANDIDATE_SPACE_ARGUMENTS = 'candidate_space_arguments'
+METHOD_SCHEME = 'method_scheme'
+PREVIOUS_METHODS = 'previous_methods'
+NEXT_METHOD = 'next_method'
+PREDECESSOR_MODEL = 'predecessor_model'
+
 MODEL = 'model'
+MODELS = 'models'
 
 # Parameters can be fixed to a value, or estimated if indicated with the string
 # `ESTIMATE`.
@@ -111,14 +118,20 @@
 
 class Method(str, Enum):
     """String literals for model selection methods."""
+
     BACKWARD = 'backward'
     BIDIRECTIONAL = 'bidirectional'
     BRUTE_FORCE = 'brute_force'
+    FAMOS = 'famos'
     FORWARD = 'forward'
+    FORWARD_AND_BACKWARD = 'forward_and_backward'
     LATERAL = 'lateral'
+    MOST_DISTANT = 'most_distant'
+
 
 class Criterion(str, Enum):
     """String literals for model selection criteria."""
+
     AIC = 'AIC'
     AICC = 'AICc'
     BIC = 'BIC'
@@ -132,18 +145,23 @@ class Criterion(str, Enum):
     Method.BACKWARD,
     Method.BIDIRECTIONAL,
     Method.FORWARD,
+    Method.FORWARD_AND_BACKWARD,
     Method.LATERAL,
 ]
 # Methods that require an initial model.
 INITIAL_MODEL_METHODS = [
     Method.BACKWARD,
+    Method.BIDIRECTIONAL,
     Method.FORWARD,
+    Method.FORWARD_AND_BACKWARD,
     Method.LATERAL,
 ]
 
 # Virtual initial models can be used to initialize some initial model methods.
 VIRTUAL_INITIAL_MODEL = 'virtual_initial_model'
 VIRTUAL_INITIAL_MODEL_METHODS = [
+    Method.BACKWARD,
+    Method.BIDIRECTIONAL,
     Method.FORWARD,
-    Method.BACKWARD
+    Method.FORWARD_AND_BACKWARD,
 ]
diff --git a/petab_select/criteria.py b/petab_select/criteria.py
index 750540a3..ffe6fa84 100644
--- a/petab_select/criteria.py
+++ b/petab_select/criteria.py
@@ -3,26 +3,19 @@
 import math
 
 import petab
-from petab.C import (
-    OBJECTIVE,
-    OBJECTIVE_PRIOR_PARAMETERS,
-    OBJECTIVE_PRIOR_TYPE,
-)
+from petab.C import OBJECTIVE, OBJECTIVE_PRIOR_PARAMETERS, OBJECTIVE_PRIOR_TYPE
 
-from .constants import (
-    #LH,
-    #LLH,
-    #NLLH,
-    Criterion,
+from .constants import (  # LH,; LLH,; NLLH,
     PETAB_PROBLEM,
     TYPE_CRITERION,
+    Criterion,
 )
-#from .model import Model
 
+# from .model import Model
 
 
 # use as attribute e.g. `Model.criterion_computer`?
-class CriterionComputer():
+class CriterionComputer:
     """Compute various criterion."""
 
     def __init__(
@@ -30,9 +23,24 @@ def __init__(
         model: 'petab_select.Model',
     ):
         self.model = model
+        self._petab_problem = None
+
+    @property
+    def petab_problem(self) -> petab.Problem:
+        """The PEtab problem that corresponds to the model.
+
+        Implemented as a property such that the `petab.Problem` object
+        is only constructed if explicitly requested.
+
+        Improves speed of operations on models by a lot. For example, analysis of models
+        that already have criteria computed can skip loading their PEtab problem again.
+        """
         # TODO refactor, if `petab_problem` is going to be produced here anyway, store
         #      in model instance instead, for use elsewhere (e.g. pyPESTO)
-        self.petab_problem = model.to_petab()[PETAB_PROBLEM]
+        #      i.e.: this is a property of a `Model` instance, not `CriterionComputer`
+        if self._petab_problem is None:
+            self._petab_problem = self.model.to_petab()[PETAB_PROBLEM]
+        return self._petab_problem
 
     def __call__(self, criterion: Criterion) -> float:
         """Get a criterion value.
@@ -44,7 +52,7 @@ def __call__(self, criterion: Criterion) -> float:
         Returns:
             The criterion value.
         """
-        return getattr(self, 'get_'+criterion.value.lower())()
+        return getattr(self, 'get_' + criterion.value.lower())()
 
     def get_aic(self) -> float:
         """Get the Akaike information criterion."""
@@ -98,7 +106,9 @@ def get_lh(self) -> float:
         elif nllh is not None:
             return math.exp(-1 * nllh)
 
-        raise ValueError('Please supply the likelihood (LH) or a compatible transformation. Compatible transformations: log(LH), -log(LH).')  # noqa: E501
+        raise ValueError(
+            'Please supply the likelihood (LH) or a compatible transformation. Compatible transformations: log(LH), -log(LH).'
+        )
 
     def get_n_estimated(self) -> int:
         """Get the number of estimated parameters."""
@@ -113,26 +123,33 @@ def get_n_priors(self) -> int:
         df = self.petab_problem.parameter_df
 
         # At least one of the objective prior columns should be present.
-        if not (OBJECTIVE_PRIOR_TYPE in df or OBJECTIVE_PRIOR_PARAMETERS in df):
+        if not (
+            OBJECTIVE_PRIOR_TYPE in df or OBJECTIVE_PRIOR_PARAMETERS in df
+        ):
             return 0
 
         # If both objective prior columns are not present, raise an error.
-        if not (OBJECTIVE_PRIOR_TYPE in df and OBJECTIVE_PRIOR_PARAMETERS in df):
-            raise NotImplementedError('Currently expect that prior types are specified with prior parameters (no default values). Please provide an example for implementation.')  # noqa: E501
+        if not (
+            OBJECTIVE_PRIOR_TYPE in df and OBJECTIVE_PRIOR_PARAMETERS in df
+        ):
+            raise NotImplementedError(
+                'Currently expect that prior types are specified with prior parameters (no default values). Please provide an example for implementation.'
+            )
 
         # Expect that the number of non-empty values in both objective prior columns
         # are the same.
         if not (
             df[OBJECTIVE_PRIOR_TYPE].notna().sum()
-            ==
-            df[OBJECTIVE_PRIOR_PARAMETERS].notna().sum()
+            == df[OBJECTIVE_PRIOR_PARAMETERS].notna().sum()
         ):
-            raise NotImplementedError('Some objective prior values are missing.')
+            raise NotImplementedError(
+                'Some objective prior values are missing.'
+            )
 
         number_of_priors = df[OBJECTIVE_PRIOR_TYPE].notna().sum()
         return number_of_priors
 
-    #def get_criterion(self, id: str) -> TYPE_CRITERION:
+    # def get_criterion(self, id: str) -> TYPE_CRITERION:
     #    """Get a criterion value, by criterion ID.
     #    FIXME: superseded by `__call__`
 
@@ -163,7 +180,7 @@ def calculate_aic(
     Returns:
         The AIC value.
     """
-    return 2*(n_estimated + nllh)
+    return 2 * (n_estimated + nllh)
 
 
 def calculate_aicc(
@@ -188,11 +205,9 @@ def calculate_aicc(
     Returns:
         The AICc value.
     """
-    return (
-        calculate_aic(n_estimated, nllh)
-        + 2*n_estimated*(n_estimated + 1)
-        / (n_measurements + n_priors - n_estimated - 1)
-    )
+    return calculate_aic(n_estimated, nllh) + 2 * n_estimated * (
+        n_estimated + 1
+    ) / (n_measurements + n_priors - n_estimated - 1)
 
 
 def calculate_bic(
@@ -217,4 +232,4 @@ def calculate_bic(
     Returns:
         The BIC value.
     """
-    return n_estimated*math.log(n_measurements + n_priors) + 2*nllh
+    return n_estimated * math.log(n_measurements + n_priors) + 2 * nllh
diff --git a/petab_select/handlers.py b/petab_select/handlers.py
index 94407718..50daa78e 100644
--- a/petab_select/handlers.py
+++ b/petab_select/handlers.py
@@ -10,7 +10,7 @@
 # TODO exclusions handler
 
 
-class LimitHandler():
+class LimitHandler:
     """A handler for classes that have a limit.
 
     Attributes:
@@ -50,4 +50,3 @@ def get_limit(self) -> TYPE_LIMIT:
             The limit.
         """
         return self.limit
-
diff --git a/petab_select/misc.py b/petab_select/misc.py
index d10cb3e8..0c1b8f9c 100644
--- a/petab_select/misc.py
+++ b/petab_select/misc.py
@@ -1,12 +1,12 @@
 import hashlib
-#import json
+
+# import json
 from typing import Any, Dict, List, Union
 
-from .constants import (
+from .constants import (  # TYPE_PARAMETER_OPTIONS_DICT,
     ESTIMATE,
     TYPE_PARAMETER_DICT,
     TYPE_PARAMETER_OPTIONS,
-    #TYPE_PARAMETER_OPTIONS_DICT,
 )
 
 
@@ -21,7 +21,7 @@ def hashify(x: Any) -> str:
     Returns:
         The hash.
     """
-    #return int(hashlib.sha256(str(x).encode('utf-8')).hexdigest(), 16)
+    # return int(hashlib.sha256(str(x).encode('utf-8')).hexdigest(), 16)
     return hashlib.blake2b(str(x).encode('utf-8')).hexdigest()
 
 
@@ -89,10 +89,9 @@ def parameter_string_to_value(
     return float_value
 
 
-
-#def hash_dictionary(dictionary: Dict[str, Union[]]):
+# def hash_dictionary(dictionary: Dict[str, Union[]]):
 #    return hash(json.dumps(dictionary, sort_keys=True))
 #
 #
-#def hash_list(list_: List):
+# def hash_list(list_: List):
 #    return hash(json.dumps(list_, sort_keys=True))
diff --git a/petab_select/model.py b/petab_select/model.py
index 6e3665fb..9a0e509d 100644
--- a/petab_select/model.py
+++ b/petab_select/model.py
@@ -1,21 +1,17 @@
 """The `Model` class."""
 import abc
-from pathlib import Path
+import warnings
 from os.path import relpath
+from pathlib import Path
 from typing import Any, Dict, List, Optional, Tuple, Union
-import warnings
-import yaml
 
-from more_itertools import one
 import numpy as np
 import petab
-from petab.C import (
-    ESTIMATE,
-    NOMINAL_VALUE,
-)
+import yaml
+from more_itertools import one
+from petab.C import ESTIMATE, NOMINAL_VALUE
 
 from .constants import (
-    Criterion,
     CRITERIA,
     ESTIMATED_PARAMETERS,
     MODEL_HASH,
@@ -28,16 +24,16 @@
     PETAB_YAML,
     PREDECESSOR_MODEL_HASH,
     TYPE_CRITERION,
-    TYPE_PATH,
     TYPE_PARAMETER,
+    TYPE_PATH,
+    VIRTUAL_INITIAL_MODEL,
+    Criterion,
 )
-from .criteria import (
-    CriterionComputer,
-)
+from .criteria import CriterionComputer
 from .misc import (
     hash_list,
-    hash_str,
     hash_parameter_dict,
+    hash_str,
     parameter_string_to_value,
 )
 from .petab import PetabMixin
@@ -110,7 +106,7 @@ class Model(PetabMixin):
     converters_save = {
         MODEL_ID: lambda x: x,
         MODEL_SUBSPACE_ID: lambda x: x,
-        MODEL_SUBSPACE_INDICES: lambda x: x,
+        MODEL_SUBSPACE_INDICES: lambda x: [int(i) for i in x],
         MODEL_HASH: lambda x: x,
         PREDECESSOR_MODEL_HASH: lambda x: x,
         PETAB_YAML: lambda x: str(x),
@@ -119,10 +115,9 @@ class Model(PetabMixin):
         # to avoid `float` cast here. Reason for cast is because e.g. pyPESTO
         # can provide type `np.float64`, which causes issues when writing to
         # YAML.
-        #ESTIMATED_PARAMETERS: lambda x: x,
+        # ESTIMATED_PARAMETERS: lambda x: x,
         ESTIMATED_PARAMETERS: lambda x: {
-            id: float(value)
-            for id, value in x.items()
+            id: float(value) for id, value in x.items()
         },
         CRITERIA: lambda x: {
             criterion_id.value: criterion_value
@@ -147,7 +142,7 @@ class Model(PetabMixin):
         #      up being calibrated to be the same... probably should be the former.
         #      Currently, the hash is stored, hence will "persist" after calibration
         #      if the same `Model` instance is used.
-        #PETAB_YAML: lambda x: hash(x),
+        # PETAB_YAML: lambda x: hash(x),
         PETAB_YAML: hash_str,
         PARAMETERS: hash_parameter_dict,
     }
@@ -210,7 +205,7 @@ def set_criterion(self, criterion: Criterion, value: float) -> None:
                 print(
                     'Overwriting saved criterion value. '
                     f'Criterion: {criterion}. Value: {self.get_criterion(criterion)}.'
-                   )
+                )
                 breakpoint()
         self.criteria[criterion] = value
 
@@ -246,8 +241,8 @@ def get_criterion(
         """
         if criterion not in self.criteria and compute:
             self.compute_criterion(criterion=criterion)
-            #value = self.criterion_computer(criterion=id)
-            #self.set_criterion(id=id, value=value)
+            # value = self.criterion_computer(criterion=id)
+            # self.set_criterion(id=id, value=value)
 
         return self.criteria.get(criterion, None)
 
@@ -272,6 +267,7 @@ def compute_criterion(self, criterion: Criterion) -> TYPE_CRITERION:
     def from_dict(
         model_dict: Dict[str, Any],
         base_path: TYPE_PATH = None,
+        petab_problem: petab.Problem = None,
     ) -> 'Model':
         """Generate a model from a dictionary of attributes.
 
@@ -285,6 +281,11 @@ def from_dict(
                 The path that any relative paths in the model are relative to
                 (e.g. the path to the PEtab problem YAML file
                 `Model.petab_yaml` may be relative).
+            petab_problem:
+                Optionally provide the PEtab problem, to avoid loading it multiple
+                times.
+                NB: This may causes issues if multiple models write to the same PEtab
+                problem in memory.
 
         Returns:
             A model instance, initialized with the provided attributes.
@@ -292,8 +293,7 @@ def from_dict(
         unknown_attributes = set(model_dict).difference(Model.converters_load)
         if unknown_attributes:
             warnings.warn(
-                'Ignoring unknown attributes: ' +
-                ', '.join(unknown_attributes)
+                'Ignoring unknown attributes: ' + ', '.join(unknown_attributes)
             )
 
         if base_path is not None:
@@ -304,6 +304,7 @@ def from_dict(
             for attribute, value in model_dict.items()
             if attribute in Model.converters_load
         }
+        model_dict[PETAB_PROBLEM] = petab_problem
         return Model(**model_dict)
 
     @staticmethod
@@ -363,8 +364,9 @@ def to_dict(
         # TODO test
         if resolve_paths:
             if model_dict[PETAB_YAML]:
-                model_dict[PETAB_YAML] = \
-                    str(Path(model_dict[PETAB_YAML]).resolve())
+                model_dict[PETAB_YAML] = str(
+                    Path(model_dict[PETAB_YAML]).resolve()
+                )
         if paths_relative_to is not None:
             if model_dict[PETAB_YAML]:
                 model_dict[PETAB_YAML] = relpath(
@@ -388,7 +390,7 @@ def to_yaml(self, petab_yaml: TYPE_PATH, *args, **kwargs) -> None:
         # kind of fixed, as the path will be resolved in `to_dict`.
         with open(petab_yaml, 'w') as f:
             yaml.dump(self.to_dict(*args, **kwargs), f)
-        #yaml.dump(self.to_dict(), str(petab_yaml))
+        # yaml.dump(self.to_dict(), str(petab_yaml))
 
     def to_petab(
         self,
@@ -406,7 +408,7 @@ def to_petab(
             with a PEtab-compatible tool for calibration of this model. If
             `output_path` is not `None`, the second value is the path to a
             PEtab YAML file that can be used to load the PEtab problem (the
-            first value) into any PEtab-compatible tool. If 
+            first value) into any PEtab-compatible tool. If
         """
         # TODO could use `copy.deepcopy(self.petab_problem)` from PetabMixin?
         petab_problem = petab.Problem.from_yaml(str(self.petab_yaml))
@@ -417,15 +419,18 @@ def to_petab(
             # Else the parameter is to be fixed.
             else:
                 petab_problem.parameter_df.loc[parameter_id, ESTIMATE] = 0
-                petab_problem.parameter_df.loc[parameter_id, NOMINAL_VALUE] = \
-                    parameter_string_to_value(parameter_value)
-                #parameter_value
+                petab_problem.parameter_df.loc[
+                    parameter_id, NOMINAL_VALUE
+                ] = parameter_string_to_value(parameter_value)
+                # parameter_value
 
         petab_yaml = None
         if output_path is not None:
             output_path = Path(output_path)
             output_path.mkdir(exist_ok=True, parents=True)
-            petab_yaml = petab_problem.to_files_generic(prefix_path=output_path)
+            petab_yaml = petab_problem.to_files_generic(
+                prefix_path=output_path
+            )
 
         return {
             PETAB_PROBLEM: petab_problem,
@@ -447,10 +452,12 @@ def get_hash(self) -> int:
             The hash.
         """
         if self.model_hash is None:
-            self.model_hash = hash_list([
-                method(getattr(self, attribute))
-                for attribute, method in Model.hash_attributes.items()
-            ])
+            self.model_hash = hash_list(
+                [
+                    method(getattr(self, attribute))
+                    for attribute, method in Model.hash_attributes.items()
+                ]
+            )
         return self.model_hash
 
     def __hash__(self) -> None:
@@ -466,9 +473,11 @@ def __str__(self):
         parameter_ids = '\t'.join(self.parameters.keys())
         parameter_values = '\t'.join(str(v) for v in self.parameters.values())
         header = '\t'.join([MODEL_ID, PETAB_YAML, parameter_ids])
-        data = '\t'.join([self.model_id, str(self.petab_yaml), parameter_values])
-        #header = f'{MODEL_ID}\t{PETAB_YAML}\t{parameter_ids}'
-        #data = f'{self.model_id}\t{self.petab_yaml}\t{parameter_values}'
+        data = '\t'.join(
+            [self.model_id, str(self.petab_yaml), parameter_values]
+        )
+        # header = f'{MODEL_ID}\t{PETAB_YAML}\t{parameter_ids}'
+        # data = f'{self.model_id}\t{self.petab_yaml}\t{parameter_values}'
         return f'{header}\n{data}'
 
     def get_mle(self) -> Dict[str, float]:
@@ -512,7 +521,10 @@ def get_estimated_parameter_ids_all(self) -> List[str]:
         # Add all estimated parameters in the PEtab problem.
         petab_problem = petab.Problem.from_yaml(str(self.petab_yaml))
         for parameter_id in petab_problem.parameter_df.index:
-            if petab_problem.parameter_df.loc[parameter_id, ESTIMATE] == PETAB_ESTIMATE_TRUE:  # noqa: E501
+            if (
+                petab_problem.parameter_df.loc[parameter_id, ESTIMATE]
+                == PETAB_ESTIMATE_TRUE
+            ):
                 estimated_parameter_ids.append(parameter_id)
 
         # Add additional estimated parameters, and collect fixed parameters,
@@ -521,8 +533,7 @@ def get_estimated_parameter_ids_all(self) -> List[str]:
         for parameter_id, value in self.parameters.items():
             if (
                 value == ESTIMATE
-                and
-                parameter_id not in estimated_parameter_ids
+                and parameter_id not in estimated_parameter_ids
             ):
                 estimated_parameter_ids.append(parameter_id)
             elif value != ESTIMATE:
@@ -562,7 +573,7 @@ def get_parameter_values(
             self.parameters.get(
                 parameter_id,
                 # Default to PEtab problem.
-                self.petab_parameters[parameter_id]
+                self.petab_parameters[parameter_id],
             )
             for parameter_id in parameter_ids
         ]
@@ -577,7 +588,7 @@ def default_compare(
     """Compare two calibrated models by their criterion values.
 
     It is assumed that the model `model0` provides a value for the criterion
-    `criterion`.
+    `criterion`, or is the `VIRTUAL_INITIAL_MODEL`.
 
     Args:
         model0:
@@ -595,10 +606,16 @@ def default_compare(
         `False`.
     """
     if not model1.has_criterion(criterion):
-        warnings.warn(f'Model "{model1.model_id}" does not provide a value for the criterion "{criterion}".')  # noqa: E501
+        warnings.warn(
+            f'Model "{model1.model_id}" does not provide a value for the criterion "{criterion}".'
+        )
         return False
+    if model0 == VIRTUAL_INITIAL_MODEL or model0 is None:
+        return True
     if criterion_threshold < 0:
-        warnings.warn('The provided criterion threshold is negative. The absolute value will be used instead.')  # noqa: E501
+        warnings.warn(
+            'The provided criterion threshold is negative. The absolute value will be used instead.'
+        )
         criterion_threshold = abs(criterion_threshold)
     if criterion in [
         Criterion.AIC,
@@ -608,8 +625,7 @@ def default_compare(
     ]:
         return (
             model1.get_criterion(criterion)
-            <
-            model0.get_criterion(criterion) - criterion_threshold
+            < model0.get_criterion(criterion) - criterion_threshold
         )
     elif criterion in [
         Criterion.LH,
@@ -617,21 +633,28 @@ def default_compare(
     ]:
         return (
             model1.get_criterion(criterion)
-            >
-            model0.get_criterion(criterion) + criterion_threshold
+            > model0.get_criterion(criterion) + criterion_threshold
         )
     else:
-        raise NotImplementedError(
-            f'Unknown criterion: {criterion}.'
-        )
+        raise NotImplementedError(f'Unknown criterion: {criterion}.')
 
 
-def models_from_yaml_list(model_list_yaml: TYPE_PATH) -> List[Model]:
+def models_from_yaml_list(
+    model_list_yaml: TYPE_PATH,
+    petab_problem: petab.Problem = None,
+    allow_single_model: bool = True,
+) -> List[Model]:
     """Generate a model from a PEtab Select list of model YAML file.
 
     Args:
         model_list_yaml:
             The path to the PEtab Select list of model YAML file.
+        petab_problem:
+            See `Model.from_dict`.
+        allow_single_model:
+            Given a YAML file that contains a single model directly (not in
+            a 1-element list), if `True` then the single model will be read in,
+            else an error will be raised.
 
     Returns:
         A list of model instances, initialized with the provided
@@ -641,8 +664,24 @@ def models_from_yaml_list(model_list_yaml: TYPE_PATH) -> List[Model]:
         model_dict_list = yaml.safe_load(f)
     if model_dict_list is None:
         return []
+
+    if not isinstance(model_dict_list, list):
+        if allow_single_model:
+            return [
+                Model.from_dict(
+                    model_dict_list,
+                    base_path=Path(model_list_yaml).parent,
+                    petab_problem=petab_problem,
+                )
+            ]
+        raise ValueError('The YAML file does not contain a list of models.')
+
     return [
-        Model.from_dict(model_dict, base_path=Path(model_list_yaml).parent)
+        Model.from_dict(
+            model_dict,
+            base_path=Path(model_list_yaml).parent,
+            petab_problem=petab_problem,
+        )
         for model_dict in model_dict_list
     ]
 
@@ -656,8 +695,7 @@ def models_to_yaml_list(
     if relative_paths:
         paths_relative_to = Path(output_yaml).parent
     model_dicts = [
-        model.to_dict(paths_relative_to=paths_relative_to)
-        for model in models
+        model.to_dict(paths_relative_to=paths_relative_to) for model in models
     ]
     model_dicts = None if not model_dicts else model_dicts
     with open(output_yaml, 'w') as f:
diff --git a/petab_select/model_space.py b/petab_select/model_space.py
index 0b898303..470dbdc6 100644
--- a/petab_select/model_space.py
+++ b/petab_select/model_space.py
@@ -1,6 +1,7 @@
 """The `ModelSpace` class and related methods."""
 import abc
 import itertools
+import logging
 from pathlib import Path
 from tempfile import NamedTemporaryFile
 from typing import (
@@ -14,10 +15,11 @@
     get_args,
 )
 
-from more_itertools import nth
 import numpy as np
 import pandas as pd
+from more_itertools import nth
 
+from .candidate_space import CandidateSpace
 from .constants import (
     ESTIMATE,
     HEADER_ROW,
@@ -31,7 +33,6 @@
     PETAB_YAML_COLUMN,
     TYPE_PATH,
 )
-from .candidate_space import CandidateSpace
 from .model import Model
 from .model_subspace import ModelSubspace
 
@@ -129,7 +130,7 @@ def line2row(
     return metadata + parameters
 
 
-class ModelSpace():
+class ModelSpace:
     """A model space, as a collection of model subspaces.
 
     Attributes:
@@ -138,6 +139,7 @@ class ModelSpace():
         exclusions:
             Hashes of models that are excluded from the model space.
     """
+
     def __init__(
         self,
         model_subspaces: List[ModelSubspace],
@@ -161,18 +163,43 @@ def from_files(
             The corresponding model space.
         """
         # TODO validate input?
-        model_space_dfs = [get_model_space_df(filename) for filename in filenames]
+        model_space_dfs = [
+            get_model_space_df(filename) for filename in filenames
+        ]
         model_subspaces = []
-        for model_space_df, model_space_filename in zip(model_space_dfs, filenames):
+        for model_space_df, model_space_filename in zip(
+            model_space_dfs, filenames
+        ):
             for model_subspace_id, definition in model_space_df.iterrows():
-                model_subspaces.append(ModelSubspace.from_definition(
+                model_subspaces.append(
+                    ModelSubspace.from_definition(
+                        model_subspace_id=model_subspace_id,
+                        definition=definition,
+                        parent_path=Path(model_space_filename).parent,
+                    )
+                )
+        model_space = ModelSpace(model_subspaces=model_subspaces)
+        return model_space
+
+    @staticmethod
+    def from_df(
+        df: pd.DataFrame,
+        parent_path: TYPE_PATH = None,
+    ):
+        model_subspaces = []
+        for model_subspace_id, definition in df.iterrows():
+            model_subspaces.append(
+                ModelSubspace.from_definition(
                     model_subspace_id=model_subspace_id,
                     definition=definition,
-                    parent_path=Path(model_space_filename).parent
-                ))
+                    parent_path=parent_path,
+                )
+            )
         model_space = ModelSpace(model_subspaces=model_subspaces)
         return model_space
 
+    # TODO: `to_df` / `to_file`
+
     def search(
         self,
         candidate_space: CandidateSpace,
@@ -194,21 +221,38 @@ def search(
             exclude:
                 Whether to exclude the new candidates from the model subspaces.
         """
-        # TODO change dict to list of subspaces. Each subspace should manage its own
-        #      ID
-        for model_subspace in self.model_subspaces.values():
-            model_subspace.search(candidate_space=candidate_space, limit=limit)
-            if len(candidate_space.models) == limit:
-                break
-            elif len(candidate_space.models) > limit:
-                raise ValueError(
-                    'An unknown error has occurred. Too many models were '
-                    f'generated. Requested limit: {limit}. Number of '
-                    f'generated models: {len(candidate_space.models)}.'
+        if candidate_space.limit.reached():
+            warnings.warn(
+                'The candidate space has already reached its limit of accepted models.',
+                RuntimeWarning,
+            )
+            return candidate_space.models
+
+        @candidate_space.wrap_search_subspaces
+        def search_subspaces(only_one_subspace: bool = False):
+            # TODO change dict to list of subspaces. Each subspace should manage its own
+            #      ID
+            if only_one_subspace and len(self.model_subspaces) > 1:
+                logging.warning(
+                    f'There is more than one model subspace. This can lead to problems for candidate space {candidate_space}, especially if they have different PEtab YAML files.'
                 )
+            for model_subspace in self.model_subspaces.values():
+                model_subspace.search(
+                    candidate_space=candidate_space, limit=limit
+                )
+                if len(candidate_space.models) == limit:
+                    break
+                elif len(candidate_space.models) > limit:
+                    raise ValueError(
+                        'An unknown error has occurred. Too many models were '
+                        f'generated. Requested limit: {limit}. Number of '
+                        f'generated models: {len(candidate_space.models)}.'
+                    )
+
+        search_subspaces()
 
         ## FIXME implement source_path.. somewhere
-        #if self.source_path is not None:
+        # if self.source_path is not None:
         #    for model in candidate_space.models:
         #        # TODO do this change elsewhere instead?
         #        # e.g. model subspace
@@ -236,7 +280,7 @@ def exclude_models(self, models: Iterable[Model]):
         # and space level.
         for model_subspace in self.model_subspaces.values():
             model_subspace.exclude_models(models)
-            #model_subspace.reset_exclusions()
+            # model_subspace.reset_exclusions()
 
     def exclude_model_hashes(self, model_hashes: Iterable[str]):
         # FIXME add Exclusions Mixin (or object) to handle exclusions on the subspace
@@ -254,10 +298,11 @@ def reset_exclusions(
 
 
 def get_model_space_df(df: Union[TYPE_PATH, pd.DataFrame]) -> pd.DataFrame:
-    #model_space_df = pd.read_csv(filename, sep='\t', index_col=MODEL_SUBSPACE_ID)  # FIXME
+    # model_space_df = pd.read_csv(filename, sep='\t', index_col=MODEL_SUBSPACE_ID)  # FIXME
     if isinstance(df, get_args(TYPE_PATH)):
         df = pd.read_csv(df, sep='\t')
-    df.set_index([MODEL_SUBSPACE_ID], inplace=True)
+    if df.index.name != MODEL_SUBSPACE_ID:
+        df.set_index([MODEL_SUBSPACE_ID], inplace=True)
     return df
 
 
@@ -265,9 +310,9 @@ def write_model_space_df(df: pd.DataFrame, filename: TYPE_PATH) -> None:
     df.to_csv(filename, sep='\t', index=True)
 
 
-#def get_model_space(
+# def get_model_space(
 #    filename: TYPE_PATH,
-#) -> List[ModelSubspace]:
+# ) -> List[ModelSubspace]:
 #    model_space_df = get_model_space_df(filename)
 #    model_subspaces = []
 #    for definition in model_space_df.iterrows():
diff --git a/petab_select/model_subspace.py b/petab_select/model_subspace.py
index e0f6361e..95efdedd 100644
--- a/petab_select/model_subspace.py
+++ b/petab_select/model_subspace.py
@@ -1,36 +1,34 @@
-from itertools import chain, product
 import math
+import warnings
+from itertools import chain, product
 from pathlib import Path
 from typing import Any, Dict, Iterable, Iterator, List, Optional, Union
-import warnings
 
-from more_itertools import one, powerset
 import numpy as np
 import pandas as pd
 import petab
-from petab.C import (
-    NOMINAL_VALUE,
-)
+from more_itertools import one, powerset
+from petab.C import NOMINAL_VALUE
 
+from .candidate_space import CandidateSpace
 from .constants import (
-    STEPWISE_METHODS,
-    Method,
     CODE_DELIMITER,
     ESTIMATE,
-    MODEL_SUBSPACE_ID,
     MODEL_SPACE_FILE_NON_PARAMETER_COLUMNS,
+    MODEL_SUBSPACE_ID,
     PARAMETER_VALUE_DELIMITER,
-    PETAB_YAML,
     PETAB_ESTIMATE_FALSE,
     PETAB_ESTIMATE_TRUE,
+    PETAB_YAML,
+    STEPWISE_METHODS,
     TYPE_PARAMETER,
-    TYPE_PARAMETER_OPTIONS,
     TYPE_PARAMETER_DICT,
+    TYPE_PARAMETER_OPTIONS,
     TYPE_PARAMETER_OPTIONS_DICT,
     TYPE_PATH,
     VIRTUAL_INITIAL_MODEL,
+    Method,
 )
-from .candidate_space import CandidateSpace
 from .misc import parameter_string_to_value
 from .model import Model
 from .petab import PetabMixin
@@ -55,6 +53,7 @@ class ModelSubspace(PetabMixin):
             Hashes of models that have been previously submitted to a candidate space
             for consideration (`CandidateSpace.consider`).
     """
+
     def __init__(
         self,
         model_subspace_id: str,
@@ -90,13 +89,9 @@ def check_compatibility_stepwise_method(
         """
         if candidate_space.method not in STEPWISE_METHODS:
             return True
-        if (
-            candidate_space.predecessor_model != VIRTUAL_INITIAL_MODEL
-            and
-            (
-                str(candidate_space.predecessor_model.petab_yaml.resolve())
-                != str(self.petab_yaml.resolve())
-            )
+        if candidate_space.predecessor_model != VIRTUAL_INITIAL_MODEL and (
+            str(candidate_space.predecessor_model.petab_yaml.resolve())
+            != str(self.petab_yaml.resolve())
         ):
             warnings.warn(
                 'The supplied candidate space is initialized with a model '
@@ -166,16 +161,15 @@ def get_models(self, estimated_parameters: List[str]) -> Iterator[Model]:
         ]
         # Generate models
         for fixed_parameter_values in product(*fixed_options):
-            fixed_parameters = dict(zip(
-                fixed_parameter_ids,
-                fixed_parameter_values,
-            ))
+            fixed_parameters = dict(
+                zip(
+                    fixed_parameter_ids,
+                    fixed_parameter_values,
+                )
+            )
             parameters = {
                 **fixed_parameters,
-                **{
-                    id: ESTIMATE
-                    for id in estimated_parameters
-                },
+                **{id: ESTIMATE for id in estimated_parameters},
             }
             model = self.parameters_to_model(parameters)
             # Skip models that are excluded.
@@ -227,13 +221,10 @@ def continue_searching(
                 return False
             return True
 
-        if not self.check_compatibility_stepwise_method(
-            candidate_space
-        ):
+        if not self.check_compatibility_stepwise_method(candidate_space):
             return
 
         if candidate_space.limit.reached():
-            warnings.warn('The candidate space has already reached its limit of accepted models.', RuntimeWarning)
             return
 
         # Compute parameter sets that are useful for finding minimal forward or backward
@@ -247,29 +238,44 @@ def continue_searching(
             else:
                 # Should already be handled elsewhere (e.g.
                 # `self.check_compatibility_stepwise_method`).
-                raise NotImplementedError(f'The default parameter set for a candidate space with the virtual initial model and method {candidate_space.method} is not implemented. Please report if this is desired.')
+                raise NotImplementedError(
+                    f'The default parameter set for a candidate space with the virtual initial model and method {candidate_space.method} is not implemented. Please report if this is desired.'
+                )
         else:
             old_estimated_all = set()
             if isinstance(candidate_space.predecessor_model, Model):
-                old_estimated_all = \
+                old_estimated_all = (
                     candidate_space.predecessor_model.get_estimated_parameter_ids_all()
+                )
 
         # Parameters that are fixed in the candidate space
         # predecessor model but are necessarily estimated in this subspace.
-        new_must_estimate_all = set(self.must_estimate_all).difference(old_estimated_all)  # noqa: E501
-        new_can_fix_all = set(old_estimated_all).difference(self.must_estimate_all)
-        new_must_fix_all = set(old_estimated_all).difference(self.can_estimate_all)
-        new_can_estimate_all = set(self.can_estimate_all).difference(old_estimated_all)
+        new_must_estimate_all = set(self.must_estimate_all).difference(
+            old_estimated_all
+        )
+        new_can_fix_all = set(old_estimated_all).difference(
+            self.must_estimate_all
+        )
+        new_must_fix_all = set(old_estimated_all).difference(
+            self.can_estimate_all
+        )
+        new_can_estimate_all = set(self.can_estimate_all).difference(
+            old_estimated_all
+        )
 
         # Parameters related to minimal changes compared to the predecessor model.
         old_estimated = set(old_estimated_all).intersection(self.can_estimate)
-        new_must_estimate = set(new_must_estimate_all).intersection(self.parameters)
+        new_must_estimate = set(new_must_estimate_all).intersection(
+            self.parameters
+        )
         # TODO remove this block...
         if not (
             set(self.must_estimate).difference(old_estimated)
             == new_must_estimate
         ):
-            raise ValueError('Unexpected error (sets that should be equal are not).')
+            raise ValueError(
+                'Unexpected error (sets that should be equal are not).'
+            )
         new_can_estimate_optional = (
             set(self.can_estimate)
             .difference(self.must_estimate)
@@ -285,8 +291,8 @@ def continue_searching(
             # There are no parameters that could become estimated in this subspace, so
             # there are no valid "forward" moves.
             if (
-                not new_can_estimate_all and
-                candidate_space.predecessor_model != VIRTUAL_INITIAL_MODEL
+                not new_can_estimate_all
+                and candidate_space.predecessor_model != VIRTUAL_INITIAL_MODEL
             ):
                 return
             # There are estimated parameters in the predecessor model that
@@ -296,26 +302,24 @@ def continue_searching(
             # Smallest possible "forward" moves involve necessarily estimated
             # parameters.
             if (
-                new_must_estimate_all or
-                candidate_space.predecessor_model == VIRTUAL_INITIAL_MODEL
+                new_must_estimate_all
+                or candidate_space.predecessor_model == VIRTUAL_INITIAL_MODEL
             ):
-                #raise NotImplementedError('untested')
                 # Consider minimal models that have all necessary parameters
                 # estimated. As this subspace necessarily has more estimated
                 # parameters than the predecessor model, all parameters that
                 # can be fixed, will be fixed.
                 estimated_parameters = {
                     parameter_id: ESTIMATE
-                    for parameter_id in
-                    set(self.parameters).difference(self.can_fix)
+                    for parameter_id in set(self.parameters).difference(
+                        self.can_fix
+                    )
                 }
                 models = self.get_models(
                     estimated_parameters=estimated_parameters,
                 )
+                previous_number_of_candidates = len(candidate_space.models)
                 for model in models:
-                    #candidate_space.consider(model)
-                    #if not continue_searching():
-                    #    return
                     continue_sending = self.send_model_to_candidate_space(
                         model=model,
                         candidate_space=candidate_space,
@@ -327,54 +331,66 @@ def continue_searching(
                 # be worse than the current set of models in the candidate
                 # space, if suitable candidates models have already been
                 # identified.
-                if candidate_space.models:
+                if len(candidate_space.models) > previous_number_of_candidates:
                     return
 
             # Keep track of the number of additional parameters estimated.
             # Stop considering models once all parameter sets with the minimal number of
             # extra estimated parameters are considered.
             n_estimated_extra = np.inf
+            previous_number_of_candidates = len(candidate_space.models)
             # The powerset should be in ascending order by number of elements.
             for parameter_set in powerset(new_can_estimate_optional):
-                # The case of a "minimal" model in the subspace being a valid candidate
-                # in this case should have been handled above already with
-                # `new_must_estimate_all`
-                if not parameter_set:
-                    continue
-                # Stop considering models as candidates once all models with a minimal
-                # increase in the number of extra estimated parameters are considered.
-                if len(parameter_set) > n_estimated_extra:
-                    break
-                n_estimated_extra = len(parameter_set)
-                estimated_parameters = (
-                    set(old_estimated)
-                    .union(new_must_estimate)
-                    .union(parameter_set)
-                )
-                models = self.get_models(
-                    estimated_parameters=list(estimated_parameters),
-                )
-                for model in models:
-                    continue_sending = self.send_model_to_candidate_space(
-                        model=model,
-                        candidate_space=candidate_space,
+                try:
+                    # The case of a "minimal" model in the subspace being a valid candidate
+                    # in this case should have been handled above already with
+                    # `new_must_estimate_all`
+                    if not parameter_set:
+                        continue
+                    # If a model has been accepted by the candidate space, only
+                    # consider models of the same size (same minimal increase in the number
+                    # of extra estimated parameters), then stop.
+                    if len(parameter_set) > n_estimated_extra:
+                        break
+                    estimated_parameters = (
+                        set(old_estimated)
+                        .union(new_must_estimate)
+                        .union(parameter_set)
                     )
-                    if not continue_searching(continue_sending):
-                        return
+                    models = self.get_models(
+                        estimated_parameters=list(estimated_parameters),
+                    )
+                    for model in models:
+                        continue_sending = self.send_model_to_candidate_space(
+                            model=model,
+                            candidate_space=candidate_space,
+                        )
+                        if not continue_searching(continue_sending):
+                            return
+                    # If model accepted set the maximal number of extra parameters to
+                    # current number of extra parameters
+                    if (
+                        len(candidate_space.models)
+                        > previous_number_of_candidates
+                    ):
+                        n_estimated_extra = len(parameter_set)
+                except StopIteration:
+                    break
 
         elif candidate_space.method == Method.BACKWARD:
             # There are no parameters that could become fixed in this subspace, so there
             # are no valid "backward" moves.
             if (
-                not new_can_fix_all and
-                candidate_space.predecessor_model != VIRTUAL_INITIAL_MODEL
+                not new_can_fix_all
+                and candidate_space.predecessor_model != VIRTUAL_INITIAL_MODEL
             ):
                 return
             # There are fixed parameters in the predecessor model that must be estimated
             # in this subspace, so there are no valid "backward" moves.
             if new_must_estimate_all:
                 return
-            # Smallest possible "backward" moves involve necessarily fixed parameters.
+            # Smallest possible "backward" moves involve necessarily fixed
+            # parameters.
             if (
                 new_must_fix_all
                 or candidate_space.predecessor_model == VIRTUAL_INITIAL_MODEL
@@ -385,12 +401,14 @@ def continue_searching(
                 # estimated.
                 estimated_parameters = {
                     parameter_id: ESTIMATE
-                    for parameter_id in
-                    set(self.parameters).difference(self.must_fix)
+                    for parameter_id in set(self.parameters).difference(
+                        self.must_fix
+                    )
                 }
                 models = self.get_models(
                     estimated_parameters=estimated_parameters,
                 )
+                previous_number_of_candidates = len(candidate_space.models)
                 for model in models:
                     continue_sending = self.send_model_to_candidate_space(
                         model=model,
@@ -403,40 +421,51 @@ def continue_searching(
                 # be worse than the current set of models in the candidate
                 # space, if suitable candidates models have already been
                 # identified.
-                if candidate_space.models:
+                if len(candidate_space.models) > previous_number_of_candidates:
                     return
 
             # Keep track of the number of new fixed parameters.
             # Stop considering models once all parameter sets with the minimal number of
             # new fixed parameters are considered.
             n_new_fixed = np.inf
+            previous_number_of_candidates = len(candidate_space.models)
             # The powerset should be in ascending order by number of elements.
             for parameter_set in powerset(new_can_fix_optional):
-                # The case of a "minimal" model in the subspace being a valid candidate
-                # in this case should have been handled above already with
-                # `new_must_estimate_all`
-                if not parameter_set:
-                    continue
-                # Stop considering models as candidates once all models with a minimal
-                # increase in the number of new fixed parameters are considered.
-                if len(parameter_set) > n_new_fixed:
-                    break
-                n_new_fixed = len(parameter_set)
-                estimated_parameters = (
-                    set(old_estimated)
-                    .union(new_must_estimate)
-                    .difference(parameter_set)
-                )
-                models = self.get_models(
-                    estimated_parameters=list(estimated_parameters),
-                )
-                for model in models:
-                    continue_sending = self.send_model_to_candidate_space(
-                        model=model,
-                        candidate_space=candidate_space,
+                try:
+                    # The case of a "minimal" model in the subspace being a valid candidate
+                    # in this case should have been handled above already with
+                    # `new_must_estimate_all`
+                    if not parameter_set:
+                        continue
+                    # If a model has been accepted by the candidate space, only
+                    # consider models of the same size (same minimal increase
+                    # in the number of new fixed parameters), then stop.
+                    if len(parameter_set) > n_new_fixed:
+                        break
+                    estimated_parameters = (
+                        set(old_estimated)
+                        .union(new_must_estimate)
+                        .difference(parameter_set)
                     )
-                    if not continue_searching(continue_sending):
-                        return
+                    models = self.get_models(
+                        estimated_parameters=list(estimated_parameters),
+                    )
+                    for model in models:
+                        continue_sending = self.send_model_to_candidate_space(
+                            model=model,
+                            candidate_space=candidate_space,
+                        )
+                        if not continue_searching(continue_sending):
+                            return
+                    # If model accepted set the number of new fixed parameters to
+                    # current number of new fixed parameters
+                    if (
+                        len(candidate_space.models)
+                        > previous_number_of_candidates
+                    ):
+                        n_new_fixed = len(parameter_set)
+                except StopIteration:
+                    break
 
         elif candidate_space.method == Method.BRUTE_FORCE:
             # TODO remove list?
@@ -453,10 +482,96 @@ def continue_searching(
                 if not continue_searching(continue_sending):
                     return
 
+        elif candidate_space.method == Method.LATERAL:
+            # There is an equal number of new necessarily estimated and fixed
+            # parameters.
+            if len(new_must_estimate_all) != len(new_must_fix_all):
+                return
+
+            if (
+                # `and` is redundant with the "equal number" check above.
+                (new_must_estimate_all and new_must_fix_all)
+                or candidate_space.predecessor_model == VIRTUAL_INITIAL_MODEL
+            ):
+                # Consider all models that have the required estimated and
+                # fixed parameters.
+                estimated_parameters = {
+                    parameter_id: ESTIMATE
+                    for parameter_id in [*old_estimated, *new_must_estimate]
+                    if parameter_id not in new_must_fix_all
+                }
+                models = self.get_models(
+                    estimated_parameters=estimated_parameters,
+                )
+                previous_number_of_candidates = len(candidate_space.models)
+                for model in models:
+                    continue_sending = self.send_model_to_candidate_space(
+                        model=model,
+                        candidate_space=candidate_space,
+                    )
+                    if not continue_searching(continue_sending):
+                        return
+
+                # No need to consider other models, as they will necessarily
+                # be worse than the current set of models in the candidate
+                # space, if suitable candidates models have already been
+                # identified.
+                if len(candidate_space.models) > previous_number_of_candidates:
+                    return
+
+            # Keep track of the number of lateral moves performed.
+            # Stop considering models once all parameter sets with the smallest
+            # lateral move size are considered.
+            n_lateral_moves = np.inf
+            previous_number_of_candidates = len(candidate_space.models)
+            # The powerset should be in ascending order by size of lateral
+            # move.
+            for parameter_set_estimate, parameter_set_fix in product(
+                powerset(new_can_estimate_optional),
+                powerset(new_can_fix_optional),
+            ):
+                try:
+                    # At least some parameters must change.
+                    if not parameter_set_estimate or not parameter_set_fix:
+                        continue
+                    # The same number of parameters must be fixed and estimated.
+                    if len(parameter_set_estimate) != len(parameter_set_fix):
+                        continue
+                    # If a model has been accepted by the candidate space, only
+                    # consider models of the same step size (same minimal
+                    # number of steps in the lateral move), then stop.
+                    if len(parameter_set_estimate) > n_lateral_moves:
+                        break
+                    estimated_parameters = (
+                        set(old_estimated)
+                        .union(new_must_estimate)
+                        .union(parameter_set_estimate)
+                        .difference(parameter_set_fix)
+                    )
+                    models = self.get_models(
+                        estimated_parameters=list(estimated_parameters),
+                    )
+                    for model in models:
+                        continue_sending = self.send_model_to_candidate_space(
+                            model=model,
+                            candidate_space=candidate_space,
+                        )
+                        if not continue_searching(continue_sending):
+                            return
+                    # If model accepted set the number of lateral moves to
+                    # current number of lateral moves
+                    if (
+                        len(candidate_space.models)
+                        > previous_number_of_candidates
+                    ):
+                        n_lateral_moves = len(parameter_set_estimate)
+                except StopIteration:
+                    break
+
         else:
             raise NotImplementedError(
                 'The requested method is not yet implemented in the model '
-                f'subspace interface: `{method}`.'
+                f'subspace interface: `{candidate_space.method}`.'
             )
 
     def send_model_to_candidate_space(
@@ -464,7 +579,7 @@ def send_model_to_candidate_space(
         model: Model,
         candidate_space: CandidateSpace,
         exclude: Optional[bool] = False,
-        #use_exclusions: Optional[bool] = True,
+        # use_exclusions: Optional[bool] = True,
     ) -> bool:
         """Send a model to a candidate space for consideration.
 
@@ -484,7 +599,7 @@ def send_model_to_candidate_space(
         # TODO if different sources of `Model` are possible (not just
         # `ModelSubspace.indices_to_model`), then would need to manage exclusions there
         # or here.
-        #if use_exclusions and hash(model) in self.exclusions:
+        # if use_exclusions and hash(model) in self.exclusions:
         #    return True
 
         if exclude:
@@ -731,8 +846,7 @@ def can_estimate_all(self) -> List[str]:
         """All parameters than can be estimated in this subspace."""
         return [
             parameter_id
-            for parameter_id, parameter_values in
-            self.parameters_all.items()
+            for parameter_id, parameter_values in self.parameters_all.items()
             if ESTIMATE in parameter_values
         ]
 
@@ -781,7 +895,6 @@ def must_estimate_all(self) -> List[str]:
         ]
         return [*must_estimate_petab, *self.must_estimate]
 
-
     def get_estimated(
         self,
         additional_parameters: Optional[TYPE_PARAMETER_DICT] = None,
@@ -809,24 +922,22 @@ def get_estimated(
 
         old_estimated_all = {
             parameter_id
-            for parameter_id, parameter_values in
-            self.parameters_all.items()
+            for parameter_id, parameter_values in self.parameters_all.items()
             if (
                 # Predecessor model sets the parameter to be estimated
                 (
-                    candidate_space
-                    .predecessor_model
-                    .parameters
-                    .get(parameter_id, None)
+                    candidate_space.predecessor_model.parameters.get(
+                        parameter_id, None
+                    )
                     == ESTIMATE
                 )
                 or (
                     # Predecessor model takes the default PEtab parameter
-                    parameter_id not in \
-                        candidate_space.predecessor_model.parameters
+                    parameter_id
+                    not in candidate_space.predecessor_model.parameters
                     and
                     # And the default PEtab parameter is estimated
-                    # The PEtab problem of this subspace and the 
+                    # The PEtab problem of this subspace and the
                     # `candidate_space` is the same, as verified earlier with
                     # `self.check_compatibility_stepwise_method`.
                     self.petab_parameters[parameter_id] == [ESTIMATE]
@@ -861,8 +972,10 @@ def decompress_parameter_values(
     parameter_strings = list(values.split(PARAMETER_VALUE_DELIMITER))
     values_decompressed = []
     for parameter_string in parameter_strings:
-        values_decompressed.append(parameter_string_to_value(
-            parameter_string=parameter_string,
-            passthrough_estimate=True,
-        ))
+        values_decompressed.append(
+            parameter_string_to_value(
+                parameter_string=parameter_string,
+                passthrough_estimate=True,
+            )
+        )
     return values_decompressed
diff --git a/petab_select/petab.py b/petab_select/petab.py
index e35f85b6..46dfb0ac 100644
--- a/petab_select/petab.py
+++ b/petab_select/petab.py
@@ -2,23 +2,17 @@
 from typing import List, Optional
 
 import petab
-from petab.C import (
-    ESTIMATE,
-    NOMINAL_VALUE,
-)
+from petab.C import ESTIMATE, NOMINAL_VALUE
 
-from .constants import (
-    PETAB_ESTIMATE_FALSE,
-    TYPE_PATH,
-    TYPE_PARAMETER_DICT,
-)
+from .constants import PETAB_ESTIMATE_FALSE, TYPE_PARAMETER_DICT, TYPE_PATH
 
 
-class PetabMixin():
+class PetabMixin:
     """Useful things for classes that contain a PEtab problem.
 
     All attributes/methods are prefixed with `petab_`.
     """
+
     def __init__(
         self,
         petab_yaml: Optional[TYPE_PATH] = None,
@@ -47,7 +41,9 @@ def __init__(
             for parameter_id, row in self.petab_problem.parameter_df.iterrows()
         }
         if parameters_as_lists:
-            self.petab_parameters = {k: [v] for k, v in self.petab_parameters.items()}
+            self.petab_parameters = {
+                k: [v] for k, v in self.petab_parameters.items()
+            }
 
     @property
     def petab_parameter_ids_estimated(self) -> List[str]:
diff --git a/petab_select/problem.py b/petab_select/problem.py
index 8147fc38..cc7695dd 100644
--- a/petab_select/problem.py
+++ b/petab_select/problem.py
@@ -3,25 +3,22 @@
 from functools import partial
 from itertools import chain
 from pathlib import Path
-from typing import Callable, Iterable, Optional, Union
+from typing import Any, Callable, Dict, Iterable, Optional, Union
+
 import yaml
 
-from .candidate_space import (
-    CandidateSpace,
-    method_to_candidate_space_class,
-)
+from .candidate_space import CandidateSpace, method_to_candidate_space_class
 from .constants import (
+    CANDIDATE_SPACE_ARGUMENTS,
     CRITERION,
     METHOD,
     MODEL_SPACE_FILES,
+    PREDECESSOR_MODEL,
     VERSION,
     Criterion,
     Method,
 )
-from .model import (
-    Model,
-    default_compare,
-)
+from .model import Model, default_compare
 from .model_space import ModelSpace
 
 
@@ -33,9 +30,12 @@ class Problem(abc.ABC):
             The model space.
         calibrated_models:
             Calibrated models. Will be used to augment the model selection problem (e.g.
-            by excluding them from the model space).
+            by excluding them from the model space). FIXME refactor out
+        candidate_space_arguments:
+            Custom options that are used to construct the candidate space.
         compare:
-            A method that compares models by selection criterion.
+            A method that compares models by selection criterion. See
+            `petab_select.model.default_compare` for an example.
         criterion:
             The criterion used to compare models.
         method:
@@ -55,9 +55,11 @@ class Problem(abc.ABC):
                 Essentially reproducible from `Problem.method` and
                 `Problem.calibrated_models`.
     """
+
     def __init__(
         self,
         model_space: ModelSpace,
+        candidate_space_arguments: Dict[str, Any] = None,
         compare: Callable[[Model, Model], bool] = None,
         criterion: Criterion = None,
         method: str = None,
@@ -70,12 +72,14 @@ def __init__(
         self.version = version
         self.yaml_path = yaml_path
 
+        self.candidate_space_arguments = candidate_space_arguments
+        if self.candidate_space_arguments is None:
+            self.candidate_space_arguments = {}
+
         self.compare = compare
         if self.compare is None:
             self.compare = partial(default_compare, criterion=self.criterion)
 
-        self.calibrated_models = []
-
     def get_path(self, relative_path: Union[str, Path]) -> Path:
         """Get the path to a resource, from a relative path.
 
@@ -94,47 +98,29 @@ def get_path(self, relative_path: Union[str, Path]) -> Path:
             return Path(relative_path)
         return self.yaml_path.parent / relative_path
 
-    def add_calibrated_models(
+    def exclude_models(
         self,
         models: Iterable[Model],
-        exclude: bool = True,
     ) -> None:
-        """Add calibrated models to the history.
+        """Exclude models from the model space.
 
         Args:
             models:
-                The models to add to the history.
-            exclude:
-                Whether to exclude the hashes of the models from the model space.
+                The models.
         """
-        self.calibrated_models = list(chain(self.calibrated_models, models))
         self.model_space.exclude_models(models)
 
-    def set_state(
+    def exclude_model_hashes(
         self,
-        state,
+        model_hashes: Iterable[str],
     ) -> None:
-        """Set the state of the problem.
-
-        Currently, only the excluded models needs to be stored.
+        """Exclude models from the model space, by model hashes.
 
         Args:
-            state:
-                The state to restore to.
-        """
-        self.add_calibrated_models(state)
-
-    def get_state(
-        self,
-    ):
-        """Get the state of the problem.
-
-        Currently, only the excluded models needs to be stored.
-
-        Returns:
-            The current state of the problem.
+            models:
+                The model hashes.
         """
-        return self.calibrated_models
+        self.model_space.exclude_model_hashes(model_hashes)
 
     @staticmethod
     def from_yaml(
@@ -160,24 +146,38 @@ def from_yaml(
             )
 
         model_space = ModelSpace.from_files(
-            #problem_specification[MODEL_SPACE_FILES],
+            # problem_specification[MODEL_SPACE_FILES],
             [
                 # `pathlib.Path` appears to handle absolute `model_space_file` paths
                 # correctly, even if used as a relative path.
                 # TODO test
                 # This is similar to the `Problem.get_path` method.
                 yaml_path.parent / model_space_file
-                for model_space_file in problem_specification[MODEL_SPACE_FILES]
+                for model_space_file in problem_specification[
+                    MODEL_SPACE_FILES
+                ]
             ],
-            #source_path=yaml_path.parent,
+            # source_path=yaml_path.parent,
         )
 
         criterion = problem_specification.get(CRITERION, None)
         if criterion is not None:
             criterion = Criterion(criterion)
 
+        candidate_space_arguments = problem_specification.get(
+            CANDIDATE_SPACE_ARGUMENTS,
+            None,
+        )
+        if candidate_space_arguments is not None:
+            if PREDECESSOR_MODEL in candidate_space_arguments:
+                candidate_space_arguments[PREDECESSOR_MODEL] = (
+                    yaml_path.parent
+                    / candidate_space_arguments[PREDECESSOR_MODEL]
+                )
+
         return Problem(
             model_space=model_space,
+            candidate_space_arguments=candidate_space_arguments,
             criterion=criterion,
             # TODO refactor method to use enum
             method=problem_specification.get(METHOD, None),
@@ -187,8 +187,9 @@ def from_yaml(
 
     def get_best(
         self,
-        models: Optional[Iterable[Model]] = None,
+        models: Optional[Iterable[Model]],
         criterion: Optional[Union[str, None]] = None,
+        compute_criterion: bool = False,
     ) -> Model:
         """Get the best model from a collection of models.
 
@@ -196,39 +197,37 @@ def get_best(
 
         Args:
             models:
-                The best model will be taken from these models. Defaults to
-                `self.calibrated_models`.
+                The best model will be taken from these models.
             criterion:
                 The criterion by which models will be compared. Defaults to
                 `self.criterion` (e.g. as defined in the PEtab Select problem YAML
                 file).
+            compute_criterion:
+                Whether to try computing criterion values, if sufficient
+                information is available (e.g., likelihood and number of
+                parameters, to compute AIC).
 
         Returns:
             The best model.
         """
         if criterion is None:
             criterion = self.criterion
-        # TODO check if commenting this out broke behavior somewhere
-        #if models is not None:
-        #    self.add_calibrated_models(models)
-        # TODO refactor s.t. `self.calibrated_models is None` when empty.
-        if models is None:
-            if self.calibrated_models:
-                models = self.calibrated_models
-            else:
-                raise ValueError('There are no calibrated models in the problem, and no models were supplied.')  # noqa: E501
 
         best_model = None
         for model in models:
+            if compute_criterion and not model.has_criterion(criterion):
+                model.get_criterion(criterion)
             if best_model is None:
                 if model.has_criterion(criterion):
                     best_model = model
                 # TODO warn if criterion is not available?
                 continue
-            if self.compare(best_model, model):
+            if self.compare(best_model, model, criterion=criterion):
                 best_model = model
         if best_model is None:
-            raise KeyError(f'None of the supplied models have a value set for the criterion {criterion}.')  # noqa: E501
+            raise KeyError(
+                f'None of the supplied models have a value set for the criterion {criterion}.'
+            )
         return best_model
 
     def new_candidate_space(
@@ -248,5 +247,17 @@ def new_candidate_space(
         if method is None:
             method = self.method
         candidate_space_class = method_to_candidate_space_class(method)
-        candidate_space = candidate_space_class(*args, **kwargs)
+        candidate_space_arguments = (
+            candidate_space_class.read_arguments_from_yaml_dict(
+                self.candidate_space_arguments
+            )
+        )
+        candidate_space_kwargs = {
+            **candidate_space_arguments,
+            **kwargs,
+        }
+        candidate_space = candidate_space_class(
+            *args,
+            **candidate_space_kwargs,
+        )
         return candidate_space
diff --git a/petab_select/ui.py b/petab_select/ui.py
index dddc8be9..04ab8d19 100644
--- a/petab_select/ui.py
+++ b/petab_select/ui.py
@@ -1,3 +1,5 @@
+import csv
+import os.path
 from pathlib import Path
 from typing import Any, Dict, List, Optional, Tuple, Union
 
@@ -6,65 +8,129 @@
 
 from .candidate_space import CandidateSpace
 from .constants import (
+    ESTIMATE,
+    INITIAL_MODEL_METHODS,
     TYPE_PATH,
-    VIRTUAL_INITIAL_MODEL_METHODS,
+    Criterion,
+    Method,
 )
-from .model import Model
+from .model import Model, default_compare
 from .problem import Problem
 
 
 def candidates(
     problem: Problem,
     candidate_space: Optional[CandidateSpace] = None,
-    predecessor_model: Optional[Model] = None,
+    previous_predecessor_model: Optional[Model] = None,
     limit: Union[float, int] = np.inf,
     limit_sent: Union[float, int] = np.inf,
+    calibrated_models: Optional[Dict[str, Model]] = None,
+    newly_calibrated_models: Optional[Dict[str, Model]] = None,
     excluded_models: Optional[List[Model]] = None,
     excluded_model_hashes: Optional[List[str]] = None,
+    criterion: Optional[Criterion] = None,
 ) -> CandidateSpace:
     """Search the model space for candidate models.
 
+    A predecessor model is chosen from `newly_calibrated_models` if available,
+    otherwise from `calibrated_models`, and is used for applicable methods.
+
     Args:
         problem:
             A PEtab Select problem.
         candidate_space:
             The candidate space. Defaults to a new candidate space based on the method
             defined in the problem.
-        predecessor_model:
-            The predecessor model for a compatible method. Defaults to the best model in
-            `Problem.calibrated_models`, if available.
+        previous_predecessor_model:
+            The previous predecessor model for a compatible method. This is
+            used as the predecessor model for the current iteration, if a
+            better model doesn't exist in the candidate space models.
         limit:
             The maximum number of models to add to the candidate space.
         limit_sent:
             The maximum number of models sent to the candidate space (which are possibly
             rejected and excluded).
+        calibrated_models:
+            All calibrated models in the model selection.
+        newly_calibrated_models:
+            All calibrated models in the most recent iteration of model
+            selection.
         excluded_models:
             Models that will be excluded from model subspaces during the search for
             candidates.
         excluded_model_hashes:
             Hashes of models that will be excluded from model subspaces during the
             search for candidates.
+        criterion:
+            The criterion by which models will be compared. Defaults to the criterion
+            defined in the PEtab Select problem.
 
     Returns:
         The candidate space, which contains the candidate models.
     """
-    if candidate_space is None:
-        candidate_space = problem.new_candidate_space(limit=limit)
-        if problem.calibrated_models:
-            candidate_space.exclude(problem.calibrated_models)
+    # FIXME might be difficult for a CLI tool to specify a specific predecessor
+    #       model if their candidate space has models. Need a way to empty
+    #       the candidate space of models... might be difficult with pickled
+    #       candidate space objects/arguments?
     if excluded_models is None:
         excluded_models = []
     if excluded_model_hashes is None:
         excluded_model_hashes = []
+    if calibrated_models is None:
+        calibrated_models = {}
+    if newly_calibrated_models is None:
+        newly_calibrated_models = {}
+    calibrated_models.update(newly_calibrated_models)
+    if criterion is None:
+        criterion = problem.criterion
+    if candidate_space is None:
+        candidate_space = problem.new_candidate_space(limit=limit)
+    candidate_space.exclude_hashes(calibrated_models)
+
+    # Set the predecessor model to the previous predecessor model.
+    # Set the new predecessor_model from the initial model or
+    # by calling ui.best to find the best model to jump to if
+    # this is not the first step of the search.
+    predecessor_model = previous_predecessor_model
+    if newly_calibrated_models:
+        predecessor_model = problem.get_best(
+            newly_calibrated_models.values(),
+            criterion=criterion,
+        )
+        # If the new predecessor model isn't better than the previous one,
+        # keep the previous one.
+        # If FAMoS jumped this will not be useful, since the jumped-to model
+        # can be expected to be worse than the jumped-from model, in general.
+        if not default_compare(
+            model0=previous_predecessor_model,
+            model1=predecessor_model,
+            criterion=criterion,
+        ):
+            predecessor_model = previous_predecessor_model
+
+        candidate_space.update_after_calibration(
+            calibrated_models=calibrated_models,
+            newly_calibrated_models=newly_calibrated_models,
+            criterion=criterion,
+        )
+        # If candidate space not Famos then ignored.
+        # Else, in case we jumped to most distant in this iteration, go into
+        # calibration with only the model we've jumped to.
+        if (
+            candidate_space.governing_method == Method.FAMOS
+            and candidate_space.jumped_to_most_distant
+        ):
+            return candidate_space
 
     if (
         predecessor_model is None
-        and
-        candidate_space.method in VIRTUAL_INITIAL_MODEL_METHODS
-        and
-        problem.calibrated_models
+        and candidate_space.method in INITIAL_MODEL_METHODS
+        and calibrated_models
     ):
-        predecessor_model = problem.get_best()
+        predecessor_model = problem.get_best(
+            models=calibrated_models.values(),
+            criterion=criterion,
+        )
     if predecessor_model is not None:
         candidate_space.reset(predecessor_model)
 
@@ -73,9 +139,18 @@ def candidates(
     #      - perhaps less portable if model IDs are generated differently on different
     #        computers
     problem.model_space.exclude_models(models=excluded_models)
-    problem.model_space.exclude_model_hashes(model_hashes=excluded_model_hashes)
+    problem.model_space.exclude_model_hashes(
+        model_hashes=excluded_model_hashes
+    )
     problem.model_space.search(candidate_space, limit=limit_sent)
 
+    write_summary_tsv(
+        problem=problem,
+        candidate_space=candidate_space,
+        previous_predecessor_model=previous_predecessor_model,
+        predecessor_model=predecessor_model,
+    )
+
     return candidate_space
 
 
@@ -144,3 +219,69 @@ def best(
     """
     # TODO return list, when multiple models are equally "best"
     return problem.get_best(models=models, criterion=criterion)
+
+
+def write_summary_tsv(
+    problem: Problem,
+    candidate_space: Optional[CandidateSpace] = None,
+    previous_predecessor_model: Optional[Union[str, Model]] = None,
+    predecessor_model: Optional[Model] = None,
+) -> None:
+    if candidate_space.summary_tsv is None:
+        return
+
+    previous_predecessor_parameter_ids = set()
+    if isinstance(previous_predecessor_model, Model):
+        previous_predecessor_parameter_ids = set(
+            previous_predecessor_model.get_estimated_parameter_ids_all()
+        )
+
+    if predecessor_model is None:
+        predecessor_model = candidate_space.predecessor_model
+    predecessor_parameter_ids = set()
+    predecessor_criterion = None
+    if isinstance(predecessor_model, Model):
+        predecessor_parameter_ids = set(
+            predecessor_model.get_estimated_parameter_ids_all()
+        )
+        predecessor_criterion = predecessor_model.get_criterion(
+            problem.criterion
+        )
+
+    diff_parameter_ids = (
+        previous_predecessor_parameter_ids.symmetric_difference(
+            predecessor_parameter_ids
+        )
+    )
+
+    diff_candidates_parameter_ids = []
+    for candidate_model in candidate_space.models:
+        candidate_parameter_ids = set(
+            candidate_model.get_estimated_parameter_ids_all()
+        )
+        diff_candidates_parameter_ids.append(
+            candidate_parameter_ids.symmetric_difference(
+                predecessor_parameter_ids
+            )
+        )
+
+    # FIXME remove once MostDistantCandidateSpace exists...
+    method = candidate_space.method
+    if (
+        candidate_space.governing_method == Method.FAMOS
+        and candidate_space.predecessor_model.predecessor_model_hash is None
+    ):
+        with open(candidate_space.summary_tsv, 'r') as f:
+            if sum(1 for _ in f) > 1:
+                method = Method.MOST_DISTANT
+
+    candidate_space.write_summary_tsv(
+        [
+            method,
+            len(candidate_space.models),
+            diff_parameter_ids,
+            predecessor_criterion,
+            predecessor_parameter_ids,
+            diff_candidates_parameter_ids,
+        ]
+    )
diff --git a/petab_select/version.py b/petab_select/version.py
index 68d146e7..8a6f55d6 100644
--- a/petab_select/version.py
+++ b/petab_select/version.py
@@ -1,2 +1,2 @@
 """Version of the model selection extension for PEtab."""
-__version__ = '0.0.9'
+__version__ = '0.1.0'
diff --git a/pyproject.toml b/pyproject.toml
index 374b58cb..50c68ba9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,3 +4,13 @@ requires = [
     "wheel"
 ]
 build-backend = "setuptools.build_meta"
+
+[tool.black]
+line-length = 79
+target-version = ['py37', 'py38', 'py39']
+skip-string-normalization = true
+
+[tool.isort]
+profile = "black"
+line_length = 79
+multi_line_output = 3
diff --git a/requirements_dev.txt b/requirements_dev.txt
new file mode 100644
index 00000000..311bf079
--- /dev/null
+++ b/requirements_dev.txt
@@ -0,0 +1,3 @@
+tox >= 3.12.4
+pre-commit >= 2.10.1
+flake8 >= 4.0.1
diff --git a/setup.py b/setup.py
index 8611e334..2faea5e1 100644
--- a/setup.py
+++ b/setup.py
@@ -1,8 +1,8 @@
-from setuptools import setup, find_packages
-import sys
 import os
 import re
+import sys
 
+from setuptools import find_packages, setup
 
 org = 'PEtab-dev'
 repo = 'petab_select'
@@ -16,10 +16,8 @@ def read(fname):
 def absolute_links(txt):
     """Replace relative petab github links by absolute links."""
 
-    raw_base = \
-        f"(https://raw.githubusercontent.com/{org}/{repo}/main/"
-    embedded_base = \
-        f"(https://github.com/{org}/{repo}/tree/main/"
+    raw_base = f"(https://raw.githubusercontent.com/{org}/{repo}/main/"
+    embedded_base = f"(https://github.com/{org}/{repo}/tree/main/"
     # iterate over links
     for var in re.findall(r'\[.*?\]\((?!http).*?\)', txt):
         if re.match(r'.*?.(png|svg)\)', var):
@@ -32,8 +30,8 @@ def absolute_links(txt):
     return txt
 
 
-# 3.7.1 for NumPy
-minimum_python_version = '3.7.1'
+# 3.7.1 for NumPy, 3.8 for `typing.get_args`
+minimum_python_version = '3.8'
 if sys.version_info < tuple(map(int, minimum_python_version.split('.'))):
     sys.exit(f'PEtab Select requires Python >= {minimum_python_version}')
 
@@ -57,12 +55,12 @@ def absolute_links(txt):
     description='PEtab Select: an extension to PEtab for model selection.',
     long_description=absolute_links(read('README.md')),
     long_description_content_type="text/markdown",
-    #author='The PEtab Select developers',
-    #author_email='dilan.pathirana@uni-bonn.de',
+    # author='The PEtab Select developers',
+    # author_email='dilan.pathirana@uni-bonn.de',
     url=f'https://github.com/{org}/{repo}',
     packages=find_packages(exclude=['doc*', 'test*']),
     install_requires=[
-        # TODO fix versions
+        # TODO minimum versions
         'more-itertools',
         'numpy',
         'pandas',
@@ -78,17 +76,17 @@ def absolute_links(txt):
         #'seaborn',
     ],
     include_package_data=True,
-    tests_require=[
-        #'flake8',
-        'pytest',
-        #'python-libcombine',
-        'amici',
-        'fides',
-        'pypesto',
-    ],
     python_requires=f'>={minimum_python_version}',
     entry_points=ENTRY_POINTS,
     extras_require={
+        'test': [
+            'pytest >= 5.4.3',
+            'pytest-cov >= 2.10.0',
+            'amici >= 0.11.25',
+            'fides >= 0.7.5',
+            #'pypesto >= 0.2.11',
+            'pypesto @ git+https://github.com/ICB-DCM/pyPESTO.git@select_method_ui_change#egg=pypesto',
+        ],
         #'reports': ['Jinja2'],
         #'combine': ['python-libcombine>=0.2.6'],
         #'doc': [
@@ -100,6 +98,6 @@ def absolute_links(txt):
         #    'nbsphinx>=0.8.2',
         #    'm2r>=0.2.1',
         #    'ipython>=7.21.0',
-        #]
-    }
+        # ]
+    },
 )
diff --git a/tox.ini b/tox.ini
new file mode 100644
index 00000000..264971a9
--- /dev/null
+++ b/tox.ini
@@ -0,0 +1,40 @@
+[tox]
+envlist =
+  clean
+  base
+
+[testenv]
+passenv = GITHUB_ACTIONS
+
+[testenv:clean]
+skip_install = true
+allowlist_externals =
+  rm
+deps = coverage
+commands =
+  coverage erase
+description =
+  Clean up
+
+[testenv:base]
+extras = test
+commands =
+  pytest --cov=petab_select --cov-report=xml --cov-append test -s
+  coverage report
+description =
+  Test basic functionality
+
+[testenv:flake8]
+skip_install = true
+deps =
+    flake8 >= 3.8.3
+    flake8-bandit >= 2.1.2
+    flake8-bugbear >= 20.1.4
+    flake8-colors >= 0.1.6
+    flake8-comprehensions >= 3.2.3
+    flake8-print >= 3.1.4
+    flake8-docstrings >= 1.6.0
+commands =
+    flake8 petab_select test setup.py
+description =
+    Run flake8 with various plugins.


Test patch:
diff --git a/tests/.gitignore b/test/.gitignore
similarity index 100%
rename from tests/.gitignore
rename to test/.gitignore
diff --git a/test/__init__.py b/test/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/test/candidate_space/__init__.py b/test/candidate_space/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/test/candidate_space/input/famos_synthetic/petab/FAMoS_2019.xml b/test/candidate_space/input/famos_synthetic/petab/FAMoS_2019.xml
new file mode 100644
index 00000000..2dbe96b4
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/petab/FAMoS_2019.xml
@@ -0,0 +1,238 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<sbml xmlns="http://www.sbml.org/sbml/level3/version1/core" level="3" version="1">
+  <model id="FAMoS_2019" name="FAMoS_2019">
+    <listOfCompartments>
+      <compartment id="Compartment" size="1" constant="true"/>
+    </listOfCompartments>
+    <listOfSpecies>
+      <species id="A" compartment="Compartment" initialAmount="100" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+      <species id="B" compartment="Compartment" initialAmount="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+      <species id="C" compartment="Compartment" initialAmount="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+      <species id="D" compartment="Compartment" initialAmount="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+    </listOfSpecies>
+    <listOfParameters>
+      <parameter id="ro_A" name="ro_A" value="0" constant="true"/>
+      <parameter id="ro_B" name="ro_B" value="0.1" constant="true"/>
+      <parameter id="ro_C" name="ro_C" value="0.1" constant="true"/>
+      <parameter id="ro_D" name="ro_D" value="0" constant="true"/>
+      <parameter id="mu_AB" name="mu_AB" value="0.1" constant="true"/>
+      <parameter id="mu_BA" name="mu_BA" value="0" constant="true"/>
+      <parameter id="mu_AC" name="mu_AC" value="0" constant="true"/>
+      <parameter id="mu_CA" name="mu_CA" value="0" constant="true"/>
+      <parameter id="mu_AD" name="mu_AD" value="0" constant="true"/>
+      <parameter id="mu_DA" name="mu_DA" value="0" constant="true"/>
+      <parameter id="mu_BC" name="mu_BC" value="0.05" constant="true"/>
+      <parameter id="mu_CB" name="mu_CB" value="0" constant="true"/>
+      <parameter id="mu_BD" name="mu_BD" value="0.2" constant="true"/>
+      <parameter id="mu_DB" name="mu_DB" value="0" constant="true"/>
+      <parameter id="mu_CD" name="mu_CD" value="0" constant="true"/>
+      <parameter id="mu_DC" name="mu_DC" value="0" constant="true"/>
+    </listOfParameters>
+    <listOfRules>
+      <rateRule variable="A">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <minus/>
+            <apply>
+              <minus/>
+              <apply>
+                <minus/>
+                <apply>
+                  <plus/>
+                  <apply>
+                    <times/>
+                    <ci> mu_BA </ci>
+                    <ci> B </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> mu_CA </ci>
+                    <ci> C </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> mu_DA </ci>
+                    <ci> D </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> ro_A </ci>
+                    <ci> A </ci>
+                  </apply>
+                </apply>
+                <apply>
+                  <times/>
+                  <ci> mu_AB </ci>
+                  <ci> A </ci>
+                </apply>
+              </apply>
+              <apply>
+                <times/>
+                <ci> mu_AC </ci>
+                <ci> A </ci>
+              </apply>
+            </apply>
+            <apply>
+              <times/>
+              <ci> mu_AD </ci>
+              <ci> A </ci>
+            </apply>
+          </apply>
+        </math>
+      </rateRule>
+      <rateRule variable="B">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <minus/>
+            <apply>
+              <minus/>
+              <apply>
+                <minus/>
+                <apply>
+                  <plus/>
+                  <apply>
+                    <times/>
+                    <ci> mu_AB </ci>
+                    <ci> A </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> mu_CB </ci>
+                    <ci> C </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> mu_DB </ci>
+                    <ci> D </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> ro_B </ci>
+                    <ci> B </ci>
+                  </apply>
+                </apply>
+                <apply>
+                  <times/>
+                  <ci> mu_BA </ci>
+                  <ci> B </ci>
+                </apply>
+              </apply>
+              <apply>
+                <times/>
+                <ci> mu_BC </ci>
+                <ci> B </ci>
+              </apply>
+            </apply>
+            <apply>
+              <times/>
+              <ci> mu_BD </ci>
+              <ci> B </ci>
+            </apply>
+          </apply>
+        </math>
+      </rateRule>
+      <rateRule variable="C">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <minus/>
+            <apply>
+              <minus/>
+              <apply>
+                <minus/>
+                <apply>
+                  <plus/>
+                  <apply>
+                    <times/>
+                    <ci> mu_AC </ci>
+                    <ci> A </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> mu_BC </ci>
+                    <ci> B </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> mu_DC </ci>
+                    <ci> D </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> ro_C </ci>
+                    <ci> C </ci>
+                  </apply>
+                </apply>
+                <apply>
+                  <times/>
+                  <ci> mu_CA </ci>
+                  <ci> C </ci>
+                </apply>
+              </apply>
+              <apply>
+                <times/>
+                <ci> mu_CB </ci>
+                <ci> C </ci>
+              </apply>
+            </apply>
+            <apply>
+              <times/>
+              <ci> mu_CD </ci>
+              <ci> C </ci>
+            </apply>
+          </apply>
+        </math>
+      </rateRule>
+      <rateRule variable="D">
+        <math xmlns="http://www.w3.org/1998/Math/MathML">
+          <apply>
+            <minus/>
+            <apply>
+              <minus/>
+              <apply>
+                <minus/>
+                <apply>
+                  <plus/>
+                  <apply>
+                    <times/>
+                    <ci> mu_AD </ci>
+                    <ci> A </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> mu_BD </ci>
+                    <ci> B </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> mu_CD </ci>
+                    <ci> C </ci>
+                  </apply>
+                  <apply>
+                    <times/>
+                    <ci> ro_D </ci>
+                    <ci> D </ci>
+                  </apply>
+                </apply>
+                <apply>
+                  <times/>
+                  <ci> mu_DA </ci>
+                  <ci> D </ci>
+                </apply>
+              </apply>
+              <apply>
+                <times/>
+                <ci> mu_DB </ci>
+                <ci> D </ci>
+              </apply>
+            </apply>
+            <apply>
+              <times/>
+              <ci> mu_DC </ci>
+              <ci> D </ci>
+            </apply>
+          </apply>
+        </math>
+      </rateRule>
+    </listOfRules>
+  </model>
+</sbml>
diff --git a/test/candidate_space/input/famos_synthetic/petab/FAMoS_2019_problem.yaml b/test/candidate_space/input/famos_synthetic/petab/FAMoS_2019_problem.yaml
new file mode 100644
index 00000000..e423f493
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/petab/FAMoS_2019_problem.yaml
@@ -0,0 +1,11 @@
+format_version: 1
+parameter_file: parameters_FAMoS_2019.tsv
+problems:
+- condition_files:
+  - experimental_conditions_FAMoS_2019.tsv
+  measurement_files:
+  - FAMos_2019_synth_measurements_0.1noise.tsv
+  observable_files:
+  - observables_FAMoS_2019.tsv
+  sbml_files:
+  - FAMoS_2019.xml
diff --git a/test/candidate_space/input/famos_synthetic/petab/FAMos_2019_synth_measurements_0.1noise.tsv b/test/candidate_space/input/famos_synthetic/petab/FAMos_2019_synth_measurements_0.1noise.tsv
new file mode 100644
index 00000000..6a234da3
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/petab/FAMos_2019_synth_measurements_0.1noise.tsv
@@ -0,0 +1,121 @@
+observableId	preequilibrationConditionId	simulationConditionId	measurement	time
+obs_A		model1_data1	90.4052240225917	1
+obs_A		model1_data1	81.9155490462083	2
+obs_A		model1_data1	74.118638599112	3
+obs_A		model1_data1	67.076299526445	4
+obs_A		model1_data1	60.7313867042841	5
+obs_A		model1_data1	54.9758314698805	6
+obs_A		model1_data1	49.6249591089255	7
+obs_A		model1_data1	45.0169961394269	8
+obs_A		model1_data1	40.6823757392052	9
+obs_A		model1_data1	36.9752363960308	10
+obs_A		model1_data1	33.2745079749387	11
+obs_A		model1_data1	30.0742660095655	12
+obs_A		model1_data1	27.3931545498098	13
+obs_A		model1_data1	24.670937708162	14
+obs_A		model1_data1	22.1726221101129	15
+obs_A		model1_data1	20.2857836657795	16
+obs_A		model1_data1	18.3124205334985	17
+obs_A		model1_data1	16.4529450566483	18
+obs_A		model1_data1	14.9595322104523	19
+obs_A		model1_data1	13.6178352834443	20
+obs_A		model1_data1	12.2698759532401	21
+obs_A		model1_data1	10.9117793381776	22
+obs_A		model1_data1	10.0075269655117	23
+obs_A		model1_data1	9.12894000775891	24
+obs_A		model1_data1	8.17552471715755	25
+obs_A		model1_data1	7.58832637116343	26
+obs_A		model1_data1	6.81481747306769	27
+obs_A		model1_data1	6.01069041893428	28
+obs_A		model1_data1	5.4656594358222	29
+obs_A		model1_data1	5.09509789405046	30
+obs_B		model1_data1	9.0007018297756	1
+obs_B		model1_data1	15.6885551330214	2
+obs_B		model1_data1	20.5555275008022	3
+obs_B		model1_data1	24.3505896189885	4
+obs_B		model1_data1	26.8387134615526	5
+obs_B		model1_data1	28.5579108201491	6
+obs_B		model1_data1	29.3386378993961	7
+obs_B		model1_data1	29.6842414879669	8
+obs_B		model1_data1	29.6012601902716	9
+obs_B		model1_data1	29.1563989173142	10
+obs_B		model1_data1	28.0915597379842	11
+obs_B		model1_data1	27.1771823929323	12
+obs_B		model1_data1	26.1996449432414	13
+obs_B		model1_data1	24.8894795720939	14
+obs_B		model1_data1	23.7958092736297	15
+obs_B		model1_data1	22.2750402794255	16
+obs_B		model1_data1	20.9051510979545	17
+obs_B		model1_data1	19.60915015308	18
+obs_B		model1_data1	18.3055488298787	19
+obs_B		model1_data1	17.1316305109407	20
+obs_B		model1_data1	15.9130210671752	21
+obs_B		model1_data1	14.6561680021543	22
+obs_B		model1_data1	13.7379948520348	23
+obs_B		model1_data1	12.688188384814	24
+obs_B		model1_data1	11.6886236992598	25
+obs_B		model1_data1	10.9831245313639	26
+obs_B		model1_data1	10.0710098586007	27
+obs_B		model1_data1	9.27665785065283	28
+obs_B		model1_data1	8.47256283413784	29
+obs_B		model1_data1	7.71049462935505	30
+obs_C		model1_data1	0.232989094821477	1
+obs_C		model1_data1	0.993499632351299	2
+obs_C		model1_data1	2.08965948178993	3
+obs_C		model1_data1	3.45518423336059	4
+obs_C		model1_data1	4.95084252305398	5
+obs_C		model1_data1	7.06264408709022	6
+obs_C		model1_data1	9.2858902955231	7
+obs_C		model1_data1	11.7393663816763	8
+obs_C		model1_data1	14.606714600774	9
+obs_C		model1_data1	17.6742934532886	10
+obs_C		model1_data1	20.9938653132865	11
+obs_C		model1_data1	24.5937563937411	12
+obs_C		model1_data1	28.7675196186448	13
+obs_C		model1_data1	33.0775441425221	14
+obs_C		model1_data1	37.8341206927535	15
+obs_C		model1_data1	43.1532681672162	16
+obs_C		model1_data1	48.7396795490924	17
+obs_C		model1_data1	54.9712638711448	18
+obs_C		model1_data1	61.7750014382742	19
+obs_C		model1_data1	69.2063960173405	20
+obs_C		model1_data1	77.1551044736292	21
+obs_C		model1_data1	86.3278593027893	22
+obs_C		model1_data1	95.9685934407861	23
+obs_C		model1_data1	106.714805376763	24
+obs_C		model1_data1	118.817854048206	25
+obs_C		model1_data1	131.819296657749	26
+obs_C		model1_data1	146.388431918277	27
+obs_C		model1_data1	161.981657265843	28
+obs_C		model1_data1	179.473855137439	29
+obs_C		model1_data1	198.770987150463	30
+obs_D		model1_data1	0.731687120540101	1
+obs_D		model1_data1	3.42150866001713	2
+obs_D		model1_data1	7.13277164777543	3
+obs_D		model1_data1	11.3626537046561	4
+obs_D		model1_data1	16.750832125073	5
+obs_D		model1_data1	22.330317912854	6
+obs_D		model1_data1	28.0551756808811	7
+obs_D		model1_data1	33.8335955902643	8
+obs_D		model1_data1	39.9060773338183	9
+obs_D		model1_data1	45.7358426721851	10
+obs_D		model1_data1	51.3568108138562	11
+obs_D		model1_data1	57.0227486341004	12
+obs_D		model1_data1	62.3225424428748	13
+obs_D		model1_data1	67.269121150386	14
+obs_D		model1_data1	72.0226608118994	15
+obs_D		model1_data1	76.7924085667823	16
+obs_D		model1_data1	80.9736216012929	17
+obs_D		model1_data1	85.2366010964468	18
+obs_D		model1_data1	88.9322839411142	19
+obs_D		model1_data1	92.2912691178331	20
+obs_D		model1_data1	95.8866770717466	21
+obs_D		model1_data1	98.6905660097913	22
+obs_D		model1_data1	101.836946624983	23
+obs_D		model1_data1	104.3108526573	24
+obs_D		model1_data1	106.840673767566	25
+obs_D		model1_data1	108.98933023116	26
+obs_D		model1_data1	111.29006614616	27
+obs_D		model1_data1	113.111729230062	28
+obs_D		model1_data1	114.967340516731	29
+obs_D		model1_data1	116.555255981599	30
diff --git a/test/candidate_space/input/famos_synthetic/petab/experimental_conditions_FAMoS_2019.tsv b/test/candidate_space/input/famos_synthetic/petab/experimental_conditions_FAMoS_2019.tsv
new file mode 100644
index 00000000..151b28b6
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/petab/experimental_conditions_FAMoS_2019.tsv
@@ -0,0 +1,2 @@
+conditionId
+model1_data1
diff --git a/test/candidate_space/input/famos_synthetic/petab/observables_FAMoS_2019.tsv b/test/candidate_space/input/famos_synthetic/petab/observables_FAMoS_2019.tsv
new file mode 100644
index 00000000..8e28902c
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/petab/observables_FAMoS_2019.tsv
@@ -0,0 +1,5 @@
+observableId	observableFormula	noiseFormula	observableTransformation	noiseDistribution
+obs_A	A	1	lin	normal
+obs_B	B	1	lin	normal
+obs_C	C	1	lin	normal
+obs_D	D	1	lin	normal
diff --git a/test/candidate_space/input/famos_synthetic/petab/parameters_FAMoS_2019.tsv b/test/candidate_space/input/famos_synthetic/petab/parameters_FAMoS_2019.tsv
new file mode 100644
index 00000000..49efa87d
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/petab/parameters_FAMoS_2019.tsv
@@ -0,0 +1,17 @@
+parameterId	parameterScale	lowerBound	upperBound	estimate	nominalValue
+ro_A	lin	-5	5	0	0
+ro_B	lin	-5	5	1	0.1
+ro_C	lin	-5	5	1	0.1
+ro_D	lin	-5	5	0	0
+mu_AB	log10	1E-15	5	1	0.1
+mu_BA	log10	1E-15	5	0	0
+mu_AC	log10	1E-15	5	0	0
+mu_CA	log10	1E-15	5	0	0
+mu_AD	log10	1E-15	5	0	0
+mu_DA	log10	1E-15	5	0	0
+mu_BC	log10	1E-15	5	1	0.05
+mu_CB	log10	1E-15	5	0	0
+mu_BD	log10	1E-15	5	1	0.2
+mu_DB	log10	1E-15	5	0	0
+mu_CD	log10	1E-15	5	0	0
+mu_DC	log10	1E-15	5	0	0
diff --git a/test/candidate_space/input/famos_synthetic/select/FAMoS_2019_petab_select_problem.yaml b/test/candidate_space/input/famos_synthetic/select/FAMoS_2019_petab_select_problem.yaml
new file mode 100644
index 00000000..094c1dc4
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/select/FAMoS_2019_petab_select_problem.yaml
@@ -0,0 +1,59 @@
+version: beta_1
+criterion: AICc
+method: famos
+model_space_files:
+- model_space_FAMoS_2019.tsv
+candidate_space_arguments:
+  critical_parameter_sets: []
+  method_scheme:
+  - previous_methods:
+    - backward
+    - forward
+    next_method: lateral
+  - previous_methods:
+    - forward
+    - backward
+    next_method: lateral
+  - previous_methods:
+    - backward
+    - lateral
+    next_method: null
+  - previous_methods:
+    - forward
+    - lateral
+    next_method: null
+  - previous_methods:
+    - forward
+    next_method: backward
+  - previous_methods:
+    - backward
+    next_method: forward
+  - previous_methods:
+    - lateral
+    next_method: forward
+  - previous_methods:
+    - most_distant
+    next_method: lateral
+  - previous_methods: null
+    next_method: lateral
+  n_reattempts: 1
+  predecessor_model: ../test_files/predecessor_model.yaml
+  consecutive_laterals: true
+  swap_parameter_sets:
+  - - ro_A
+    - mu_BA
+    - mu_CA
+    - mu_DA
+  - - ro_B
+    - mu_AB
+    - mu_CB
+    - mu_DB
+  - - ro_C
+    - mu_AC
+    - mu_BC
+    - mu_DC
+  - - ro_D
+    - mu_AD
+    - mu_BD
+    - mu_CD
+  summary_tsv: output_famos/summary.tsv
diff --git a/test/candidate_space/input/famos_synthetic/select/model_space_FAMoS_2019.tsv b/test/candidate_space/input/famos_synthetic/select/model_space_FAMoS_2019.tsv
new file mode 100644
index 00000000..9d602c3f
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/select/model_space_FAMoS_2019.tsv
@@ -0,0 +1,2 @@
+model_subspace_id	petab_yaml	ro_A	ro_B	ro_C	ro_D	mu_AB	mu_BA	mu_AC	mu_CA	mu_AD	mu_DA	mu_BC	mu_CB	mu_BD	mu_DB	mu_CD	mu_DC
+model_subspace_1	../petab/FAMoS_2019_problem.yaml	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate	0;estimate
diff --git a/test/candidate_space/input/famos_synthetic/test_files/calibration_results.tsv b/test/candidate_space/input/famos_synthetic/test_files/calibration_results.tsv
new file mode 100644
index 00000000..4a6c0702
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/test_files/calibration_results.tsv
@@ -0,0 +1,282 @@
+model_id	AICc
+M_0000010010011011	139213.93785488105
+M_0000011010011010	147241.83250045744
+M_0000110010001011	129085.62527164364
+M_0001000011010011	20019.422755830317
+M_0001000110010011	20019.42287885885
+M_0001001011010010	58645.5159924183
+M_0001001110010010	85944.52663282724
+M_0001010000011011	561253.8821574149
+M_0001010010000111	17407.944853575715
+M_0001010010001101	11373.70265679193
+M_0001010010010011	20019.42291165012
+M_0001010010011001	19822.224390562587
+M_0001010010110010	223159.7600795708
+M_0001011000011010	119758.91035527596
+M_0001011010000110	52539.03132570987
+M_0001011010010010	85944.5264701275
+M_0001011010011000	85944.52656271843
+M_0001100011000011	8783.477653560316
+M_0001100011001001	2082.2848460706477
+M_0001100110000011	8783.477689063757
+M_0001100110001001	2122.93902388184
+M_0001110000001011	2126.584221871469
+M_0001110010000011	5989.975092399021
+M_0001110010001001	2122.9390944314687
+M_0001110010100010	66865.97525944727
+M_0001110010101000	66865.97526955212
+M_0001111010000010	81844.29973880656
+M_0001111010001000	81844.29970386633
+M_0010100000101000	13964.86529487194
+M_0010100000101100	13967.043858068218
+M_0010100000111000	4668.469165871618
+M_0011010010010010	223159.76006183933
+M_0011100000101000	2745.410237113948
+M_0011100000101001	1245.2636262725478
+M_0011110010000010	209279.1634871877
+M_0011110010001000	206824.56340461335
+M_0100100000101000	23409.18328914282
+M_0100100000101001	10168.445671589136
+M_0100101000101000	23411.36177928765
+M_0101010010000011	21480.37645074097
+M_0101010010001001	21480.37627690475
+M_0101011010000010	85944.52650071045
+M_0101100000101000	22297.504133595157
+M_0101100000101001	1828.464353816787
+M_0110000000101000	561249.48662067
+M_0110000000101100	561251.6651103725
+M_0110000000111000	561251.6651103725
+M_0110100000001000	206221.04167048345
+M_0110100000001001	2806.047402348072
+M_0110100000100000	177910.00743008128
+M_0110100000100010	33340.41218209771
+M_0110100000101000	232.113265665936
+M_0110100000101001	234.3303138148792
+M_0110100000101010	234.3156308348245
+M_0110100000101100	234.33031342028949
+M_0110100000111000	234.3268229275457
+M_0110100001101000	234.33031314833016
+M_0110100001101010	236.56805265000804
+M_0110100010100000	15756.946297002763
+M_0110100010101000	234.33031288019092
+M_0110100011101000	236.5869504485113
+M_0110100100101000	234.3273835278236
+M_0110100100101010	236.5691109524262
+M_0110100110101000	236.5840203029052
+M_0110101000001000	497.3707527496445
+M_0110101000101000	234.33031295737
+M_0110110000101000	234.32946884457883
+M_0110110000101010	236.5687357954883
+M_0110110010101000	236.58610604593983
+M_0111000000101000	561251.6651103725
+M_0111000000101001	561253.8821574149
+M_0111100000001000	206223.21384086332
+M_0111100000001001	606.9037822865472
+M_0111100000100000	177912.18591824788
+M_0111100000100001	177914.40296609284
+M_0111100000101000	234.3250077924832
+M_0111100000101001	236.54705233061625
+M_0111100000101101	238.84432917693465
+M_0111100010101000	236.5816448113885
+M_0111100010101001	238.84434958475077
+M_0111100100101101	241.1790156115896
+M_0111100101101101	243.56100083729265
+M_0111101010101001	241.18341611694956
+M_0111110101101101	245.98625534924483
+M_1000100010001011	25488.987549211623
+M_1000110111010111	33682.94076210533
+M_1001000010001101	10059.448865593136
+M_1001000010010011	14983.46222088898
+M_1001000010011001	14983.462223649214
+M_1001001010010010	62368.876932390696
+M_1001100000001011	1876.108100244288
+M_1001100010000011	7554.901601882608
+M_1001100010001001	1876.108025959749
+M_1001100010001011	1878.3646939136331
+M_1001100010001101	1874.3052802249747
+M_1001100010011001	1878.3647690190437
+M_1001100010101000	29004.845018270666
+M_1001100010101001	1878.364725287871
+M_1001100011001001	1878.364628524238
+M_1001100110001001	1878.364779043418
+M_1001101010001000	48417.47097786767
+M_1001101010001001	1878.3645949487295
+M_1001110010001001	1878.3647795773736
+M_1001110101011101	1881.3236370601333
+M_1001110111010101	5512.646485215627
+M_1010100000101010	554.6582992812749
+M_1010100000101110	506.043663839271
+M_1010100000111010	497.6509421175728
+M_1010100010101000	882.8702477394546
+M_1010100010101100	885.1269082847127
+M_1010100010111000	601.7940759062141
+M_1011100000101001	433.3882550546375
+M_1011100000101101	360.01073998316195
+M_1011100010001000	206269.22383117743
+M_1011100010001001	572.4115972853955
+M_1011100010001011	574.7089604363986
+M_1011100010001101	574.2624322321558
+M_1011100010011001	574.7073976582066
+M_1011100010101000	405.8386721143808
+M_1011100010101001	382.5778905589666
+M_1011100011001001	574.7089621871352
+M_1011100100101101	362.3498014779889
+M_1011100101101101	364.7317910137515
+M_1011100110001001	574.7089926825686
+M_1011101010001001	301.201692408097
+M_1011101010001011	303.5407915622274
+M_1011101010001101	295.47824301127844
+M_1011101010011001	301.52029337677715
+M_1011101010101001	303.54077561640514
+M_1011101011001001	303.5407901163882
+M_1011101110001001	303.54079317899004
+M_1011110010001001	574.7089840081791
+M_1011110101101101	367.15788201970594
+M_1011111010001001	303.54077122448376
+M_1100010111010111	49813.91784706038
+M_1100100000101010	17877.442288126847
+M_1100100000101011	4770.4712415596505
+M_1100100010101000	9773.734269092403
+M_1100100010101001	4770.471310151174
+M_1100101000101010	17879.698953110583
+M_1100101010101000	9775.990866941636
+M_1100110101001111	4777.489639680295
+M_1100110111000111	30332.80096997364
+M_1100110111001101	4777.489716081316
+M_1100110111010011	30332.800929159228
+M_1100110111100110	11016.328616609726
+M_1100111111000110	41228.25894861741
+M_1101000010001001	16723.60491947643
+M_1101010101011101	397001.17155116936
+M_1101010111010101	10043.947785511777
+M_1101100000101001	1828.4669039830949
+M_1101100000101101	1830.7642054619791
+M_1101100010001001	1828.466866432956
+M_1101100010101000	9512.205257685406
+M_1101100010101001	1830.764168943262
+M_1101100100101101	1833.1032679486905
+M_1101100101101101	1835.4852513465949
+M_1101101010101001	1833.1032581612865
+M_1101110101000111	383914.7604399954
+M_1101110101001101	1835.4852672346087
+M_1101110101001111	1837.9113687541485
+M_1101110101011001	1835.4852758742516
+M_1101110101011101	1837.9113595304132
+M_1101110101101100	9531.431831583666
+M_1101110101101101	1837.9113687612264
+M_1101110111000101	3646.48444481816
+M_1101110111001101	1837.9113548340729
+M_1101110111010001	3646.484428093567
+M_1101110111100100	9470.132468078638
+M_1101111101001100	41030.44544078736
+M_1101111101001101	1837.9113562254213
+M_1101111111000100	41226.41791870745
+M_1110000000101010	396991.89656571264
+M_1110000000101110	396994.1532027124
+M_1110000000111010	396994.1532025714
+M_1110000010101000	225250.80565647891
+M_1110000010101100	7354.181781450917
+M_1110000010111000	225253.0622635607
+M_1110100000001010	206225.4249527701
+M_1110100000001011	709.8139285203328
+M_1110100000100010	32569.11777402543
+M_1110100000101000	234.32053581346983
+M_1110100000101001	236.57717236056465
+M_1110100000101010	236.5309219253127
+M_1110100000101011	238.7752196664412
+M_1110100000101101	238.87446942005604
+M_1110100000101110	238.81920116626029
+M_1110100000111010	238.8209933841467
+M_1110100001101010	238.82769297452745
+M_1110100010001000	206225.42523656492
+M_1110100010001001	1903.1635775769812
+M_1110100010100000	6175.236251468945
+M_1110100010100010	2435.7230027787205
+M_1110100010101000	236.57717283177325
+M_1110100010101001	238.87447010638908
+M_1110100010101010	238.82096515868363
+M_1110100010101100	238.8744691260788
+M_1110100010111000	238.85880783638197
+M_1110100011101000	238.87446936299148
+M_1110100100101010	238.8271393711566
+M_1110100100101101	241.2116011501129
+M_1110100101101101	243.593585788594
+M_1110100110101000	238.87253451641828
+M_1110101000001010	381.953412082261
+M_1110101000101010	238.82821959285397
+M_1110101010001000	382.0059312235692
+M_1110101010101000	238.87446999761505
+M_1110101010101001	241.21353703929177
+M_1110110000101010	238.82821945073889
+M_1110110010101000	238.8744704850313
+M_1110110101101101	246.0196818158347
+M_1110110111000110	211557.76697941415
+M_1111000000101001	396994.1532026271
+M_1111000000101101	396996.4504998744
+M_1111000010101000	220481.7872639834
+M_1111000010101001	14580.346343644816
+M_1111000100101101	396998.78956616647
+M_1111000101101101	397001.1715515276
+M_1111001010101001	14505.882802401526
+M_1111010101101101	397003.59764623776
+M_1111100000001001	606.9905832634503
+M_1111100000001101	602.4411335548343
+M_1111100000100001	177916.64191094312
+M_1111100000100010	14189.991928330788
+M_1111100000100101	177918.93920798044
+M_1111100000101000	236.5470995198401
+M_1111100000101001	238.80076560298963
+M_1111100000101010	238.8136978955608
+M_1111100000101100	238.8331287054401
+M_1111100000101101	241.1222202761463
+M_1111100001101101	243.50420511898895
+M_1111100010001000	206227.63311331163
+M_1111100010001001	568.9552843715751
+M_1111100010100000	1404.2903801434873
+M_1111100010100001	1406.5384305557384
+M_1111100010101000	238.83160653758773
+M_1111100010101001	241.12867656656744
+M_1111100100001101	604.7802175312609
+M_1111100100100101	177921.2780036785
+M_1111100100101001	241.13748677613373
+M_1111100100101100	241.1720398184164
+M_1111100100101101	243.50405343595625
+M_1111100101001101	607.1622021577575
+M_1111100101100101	177923.66001179285
+M_1111100101101001	243.5194719023592
+M_1111100101101100	243.55402510830623
+M_1111100101101101	245.9301495766389
+M_1111101000101001	241.12880866953333
+M_1111101010001001	261.3039091174704
+M_1111101010001011	263.6116954813844
+M_1111101010001101	263.66190443322904
+M_1111101010011001	263.65972246373093
+M_1111101010100001	1408.87016504161
+M_1111101010101000	241.17067087762751
+M_1111101010101001	243.49951157260497
+M_1111101010101011	245.90866663807364
+M_1111101010101101	245.91762078132183
+M_1111101010111001	245.91680600377367
+M_1111101011001001	263.6858866332959
+M_1111101011101001	245.92560809043647
+M_1111101110001001	263.6849507178238
+M_1111101110101001	245.92424808908
+M_1111110001101101	245.93030121541156
+M_1111110100101101	245.9301495839917
+M_1111110101001100	206234.6454773945
+M_1111110101001101	609.5826436734602
+M_1111110101001111	606.3868920612166
+M_1111110101011101	612.0540839150935
+M_1111110101100101	177926.08608462897
+M_1111110101101001	245.94556766090545
+M_1111110101101100	245.98012071080976
+M_1111110101101101	248.40159324081264
+M_1111110101101111	250.9062566610432
+M_1111110101111101	250.91785846088345
+M_1111110111000100	207402.40135257572
+M_1111110111001101	578.5738952510358
+M_1111110111101101	250.91739945255304
+M_1111111010001001	263.6858967007385
+M_1111111010101001	245.92560787196052
+M_1111111101001101	268.57688618118783
+M_1111111101101101	250.90966149938663
diff --git a/test/candidate_space/input/famos_synthetic/test_files/predecessor_model.yaml b/test/candidate_space/input/famos_synthetic/test_files/predecessor_model.yaml
new file mode 100644
index 00000000..f0442820
--- /dev/null
+++ b/test/candidate_space/input/famos_synthetic/test_files/predecessor_model.yaml
@@ -0,0 +1,39 @@
+criteria:
+  AIC: 30330.782621349786
+  AICc: 30332.80096997364
+  BIC: 30358.657538777607
+  NLLH: 15155.391310674893
+estimated_parameters:
+  mu_AB: 0.09706971737957297
+  mu_AD: -0.6055359156893474
+  mu_BA: 0.6989700040781575
+  mu_CA: -13.545121478780585
+  mu_CD: -13.955162965672203
+  mu_DA: -13.405909047226377
+  mu_DB: -13.402598631022197
+  mu_DC: -1.1619119214640863
+  ro_A: -1.6431508614147425
+  ro_B: 2.9912966824709097
+model_hash: null
+model_id: M_1100110111000111
+model_subspace_id: model_subspace_1
+model_subspace_indices: null
+parameters:
+  mu_AB: estimate
+  mu_AC: 0
+  mu_AD: estimate
+  mu_BA: estimate
+  mu_BC: 0
+  mu_BD: 0
+  mu_CA: estimate
+  mu_CB: 0
+  mu_CD: estimate
+  mu_DA: estimate
+  mu_DB: estimate
+  mu_DC: estimate
+  ro_A: estimate
+  ro_B: estimate
+  ro_C: 0
+  ro_D: 0
+petab_yaml: ../petab/FAMoS_2019_problem.yaml
+predecessor_model_hash: null
diff --git a/test/candidate_space/test_candidate_space.py b/test/candidate_space/test_candidate_space.py
new file mode 100644
index 00000000..9a993ef6
--- /dev/null
+++ b/test/candidate_space/test_candidate_space.py
@@ -0,0 +1,235 @@
+from pathlib import Path
+
+import pandas as pd
+import pytest
+from more_itertools import one
+
+import petab_select
+from petab_select.candidate_space import (  # BackwardCandidateSpace,; BruteForceCandidateSpace,; ForwardCandidateSpace,; ForwardAndBackwardCandidateSpace,; LateralCandidateSpace,
+    BidirectionalCandidateSpace,
+)
+from petab_select.constants import (
+    ESTIMATE,
+    MODEL_SUBSPACE_ID,
+    MODELS,
+    PARAMETER_VALUE_DELIMITER,
+    PETAB_YAML,
+    Criterion,
+)
+from petab_select.model import Model, default_compare
+from petab_select.model_space import ModelSpace, get_model_space_df
+from petab_select.model_subspace import ModelSubspace
+
+
+@pytest.fixture
+def ordered_model_parameterizations():
+    good_models_ascending = [
+        # forward
+        '00000',
+        '10000',
+        '11000',
+        '11100',
+        '11110',
+        # backward
+        '01110',
+        '01100',
+        # forward
+        '01101',
+        '01111',
+        # backward
+        '00111',
+        '00011',
+    ]
+    bad_models = [
+        '01011',
+        '11011',
+    ]
+
+    # All good models are unique
+    assert len(set(good_models_ascending)) == len(good_models_ascending)
+    # All bad models are unique
+    assert len(set(bad_models)) == len(bad_models)
+    # No models are defined as both bad and good.
+    assert not set(good_models_ascending).intersection(bad_models)
+
+    return good_models_ascending, bad_models
+
+
+@pytest.fixture
+def calibrated_model_space(ordered_model_parameterizations):
+    good_models_ascending, bad_models = ordered_model_parameterizations
+
+    # As good models are ordered ascending by "goodness", and criteria
+    # decreases for better models, the criteria decreases as the index increases.
+    good_model_criteria = {
+        model: 100 - index for index, model in enumerate(good_models_ascending)
+    }
+    # All bad models are currently set to the same "bad" criterion value.
+    bad_model_criteria = {model: 1000 for model in bad_models}
+
+    model_criteria = {
+        **good_model_criteria,
+        **bad_model_criteria,
+    }
+    return model_criteria
+
+
+@pytest.fixture
+def model_space(calibrated_model_space) -> pd.DataFrame:
+    data = {
+        "model_subspace_id": [],
+        "petab_yaml": [],
+        "k1": [],
+        "k2": [],
+        "k3": [],
+        "k4": [],
+        "k5": [],
+    }
+
+    for model in calibrated_model_space:
+        data["model_subspace_id"].append(f"model_subspace_{model}")
+        data["petab_yaml"].append(
+            Path(__file__).parent.parent.parent
+            / 'doc'
+            / 'examples'
+            / 'model_selection'
+            / 'petab_problem.yaml'
+        )
+        k1, k2, k3, k4, k5 = [
+            '0' if parameter == '0' else ESTIMATE for parameter in model
+        ]
+        data["k1"].append(k1)
+        data["k2"].append(k2)
+        data["k3"].append(k3)
+        data["k4"].append(k4)
+        data["k5"].append(k5)
+    df = pd.DataFrame(data=data)
+    df = get_model_space_df(df)
+    model_space = ModelSpace.from_df(df)
+    return model_space
+
+
+def test_bidirectional(
+    model_space, calibrated_model_space, ordered_model_parameterizations
+):
+    criterion = Criterion.AIC
+    model_id_length = one(
+        set([len(model_id) for model_id in calibrated_model_space])
+    )
+
+    candidate_space = BidirectionalCandidateSpace()
+    calibrated_models = []
+
+    # Perform bidirectional search until no more models are found.
+    search_iterations = 0
+    while True:
+        new_calibrated_models = []
+        search_iterations += 1
+
+        # Get models.
+        model_space.search(candidate_space)
+
+        # Calibrate models.
+        for model in candidate_space.models:
+            model_id = model.model_subspace_id[-model_id_length:]
+            model.set_criterion(
+                criterion=criterion, value=calibrated_model_space[model_id]
+            )
+            new_calibrated_models.append(model)
+
+        # End if no more models were found.
+        if not new_calibrated_models:
+            break
+
+        # Get next predecessor model as best of new models.
+        best_new_model = None
+        for model in new_calibrated_models:
+            if best_new_model is None:
+                best_new_model = model
+                continue
+            if default_compare(
+                model0=best_new_model, model1=model, criterion=criterion
+            ):
+                best_new_model = model
+
+        # Set next predecessor and exclusions.
+        calibrated_model_hashes = [
+            model.get_hash() for model in calibrated_models
+        ]
+        candidate_space.reset(
+            predecessor_model=best_new_model,
+            exclusions=calibrated_model_hashes,
+        )
+
+        # exclude calibrated model hashes from model space too?
+        model_space.exclude_model_hashes(model_hashes=calibrated_model_hashes)
+
+    # Check that expected models are found at each iteration.
+    (
+        good_model_parameterizations_ascending,
+        bad_model_parameterizations,
+    ) = ordered_model_parameterizations
+    search_iteration = 0
+    for history_item in candidate_space.method_history:
+        models = history_item[MODELS]
+        if not models:
+            continue
+        model_parameterizations = [
+            model.model_subspace_id[-5:] for model in models
+        ]
+
+        good_model_parameterization = good_model_parameterizations_ascending[
+            search_iteration
+        ]
+        # The expected good model was found.
+        assert good_model_parameterization in model_parameterizations
+        model_parameterizations.remove(good_model_parameterization)
+
+        if search_iteration == 0:
+            # All parameterizations have been correctly identified and removed.
+            assert not model_parameterizations
+            search_iteration += 1
+            continue
+
+        previous_model_parameterization = (
+            good_model_parameterizations_ascending[search_iteration - 1]
+        )
+
+        # The expected bad model is also found.
+        # If a bad model is the same dimension and also represents a similar stepwise move away from the previous
+        # model parameterization, it should also be in the parameterizations.
+        for bad_model_parameterization in bad_model_parameterizations:
+            # Skip if different dimensions
+            if sum(map(int, bad_model_parameterization)) != sum(
+                map(int, good_model_parameterization)
+            ):
+                continue
+            # Skip if different distances from previous model parameterization
+            if sum(
+                [
+                    a != b
+                    for a, b in zip(
+                        bad_model_parameterization,
+                        previous_model_parameterization,
+                    )
+                ]
+            ) != sum(
+                [
+                    a != b
+                    for a, b in zip(
+                        good_model_parameterization,
+                        previous_model_parameterization,
+                    )
+                ]
+            ):
+                continue
+            assert bad_model_parameterization in model_parameterizations
+            model_parameterizations.remove(bad_model_parameterization)
+
+        # All parameterizations have been correctly identified and removed.
+        assert not model_parameterizations
+        search_iteration += 1
+
+        # End test if all good models were found in the correct order.
+        if search_iteration >= len(good_model_parameterizations_ascending):
+            break
diff --git a/test/candidate_space/test_famos.py b/test/candidate_space/test_famos.py
new file mode 100644
index 00000000..b0dde6f2
--- /dev/null
+++ b/test/candidate_space/test_famos.py
@@ -0,0 +1,167 @@
+from pathlib import Path
+from typing import Tuple
+
+import numpy as np
+import pandas as pd
+import petab
+import pytest
+from more_itertools import one
+
+import petab_select
+from petab_select import ESTIMATE, FamosCandidateSpace, Method, Model
+from petab_select.constants import Criterion
+from petab_select.model import default_compare
+
+
+@pytest.fixture
+def input_path():
+    return Path(__file__).parent / "input" / "famos_synthetic"
+
+
+@pytest.fixture
+def petab_select_problem(input_path):
+    return petab_select.Problem.from_yaml(
+        input_path / "select" / "FAMoS_2019_petab_select_problem.yaml"
+    )
+
+
+@pytest.fixture
+def expected_criterion_values(input_path):
+    calibration_results = pd.read_csv(
+        input_path / "test_files" / "calibration_results.tsv",
+        sep="\t",
+    ).set_index('model_id')
+    return dict(calibration_results['AICc'])
+
+
+@pytest.fixture
+def expected_progress_list():
+    return [
+        (Method.LATERAL, set()),
+        (Method.LATERAL, {4, 15}),
+        (Method.LATERAL, {9, 13}),
+        (Method.FORWARD, set()),
+        (Method.FORWARD, {3}),
+        (Method.FORWARD, {11}),
+        (Method.BACKWARD, set()),
+        (Method.BACKWARD, {6}),
+        (Method.BACKWARD, {10}),
+        (Method.BACKWARD, {8}),
+        (Method.BACKWARD, {14}),
+        (Method.BACKWARD, {1}),
+        (Method.BACKWARD, {16}),
+        (Method.BACKWARD, {4}),
+        (Method.FORWARD, set()),
+        (Method.LATERAL, set()),
+        (Method.MOST_DISTANT, {2, 3, 4, 5, 6, 7, 9, 11, 12, 13, 15}),
+        (Method.LATERAL, {16, 7}),
+        (Method.LATERAL, {5, 12}),
+        (Method.LATERAL, {13, 15}),
+        (Method.LATERAL, {1, 6}),
+        (Method.FORWARD, set()),
+        (Method.FORWARD, {3}),
+        (Method.FORWARD, {7}),
+        (Method.FORWARD, {2}),
+        (Method.FORWARD, {11}),
+        (Method.BACKWARD, set()),
+        (Method.BACKWARD, {7}),
+        (Method.BACKWARD, {16}),
+        (Method.BACKWARD, {4}),
+        (Method.FORWARD, set()),
+        (Method.LATERAL, set()),
+        (Method.LATERAL, {9, 15}),
+        (Method.FORWARD, set()),
+        (Method.BACKWARD, set()),
+        (Method.LATERAL, set()),
+    ]
+
+
+def test_famos(
+    petab_select_problem,
+    expected_criterion_values,
+    expected_progress_list,
+):
+    def set_model_id(model):
+        model.model_id = "M_" + ''.join(
+            '1' if v == ESTIMATE else '0' for v in model.parameters.values()
+        )
+
+    def calibrate(
+        model,
+        expected_criterion_values=expected_criterion_values,
+    ) -> None:
+        model.set_criterion(
+            criterion=petab_select_problem.criterion,
+            value=expected_criterion_values[model.model_id],
+        )
+
+    def parse_summary_to_progress_list(summary_tsv: str) -> Tuple[Method, set]:
+        """Get progress information from the summary file."""
+        df_raw = pd.read_csv(summary_tsv, sep='\t')
+        df = df_raw.loc[~pd.isnull(df_raw["predecessor change"])]
+
+        parameter_list = list(
+            petab_select_problem.model_space.model_subspaces[
+                'model_subspace_1'
+            ].parameters
+        )
+
+        progress_list = []
+
+        for index, (_, row) in enumerate(df.iterrows()):
+            method = Method(row["method"])
+
+            model = {
+                1 + parameter_list.index(parameter_id)
+                for parameter_id in eval(row["current model"])
+            }
+            if index == 0:
+                model0 = model
+
+            difference = model.symmetric_difference(model0)
+            progress_list.append((method, difference))
+            model0 = model
+
+        return progress_list
+
+    progress_list = []
+    calibrated_models = {}
+    newly_calibrated_models = {}
+
+    candidate_space = petab_select_problem.new_candidate_space()
+    candidate_space.summary_tsv.unlink(missing_ok=True)
+    candidate_space._setup_summary_tsv()
+
+    with pytest.raises(StopIteration, match="No valid models found."):
+        predecessor_model = candidate_space.predecessor_model
+        while True:
+            # Save predecessor_models and find new candidates
+            if candidate_space.predecessor_model is not None:
+                previous_predecessor_model = candidate_space.predecessor_model
+            else:
+                previous_predecessor_model = one(candidate_space.models)
+            candidate_space = petab_select.ui.candidates(
+                problem=petab_select_problem,
+                candidate_space=candidate_space,
+                calibrated_models=calibrated_models,
+                newly_calibrated_models=newly_calibrated_models,
+                previous_predecessor_model=previous_predecessor_model,
+            )
+
+            # Calibrate candidate models
+            newly_calibrated_models = {}
+            for candidate_model in candidate_space.models:
+                set_model_id(candidate_model)
+                calibrate(candidate_model)
+                newly_calibrated_models[
+                    candidate_model.get_hash()
+                ] = candidate_model
+                calibrated_models.update(newly_calibrated_models)
+
+            # Stop iteration if there are no candidate models
+            if not candidate_space.models:
+                raise StopIteration("No valid models found.")
+
+    progress_list = parse_summary_to_progress_list(candidate_space.summary_tsv)
+
+    assert progress_list == expected_progress_list, progress_list
diff --git a/test/cli/__init__.py b/test/cli/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/cli/expected_output/model/conditions.tsv b/test/cli/expected_output/model/conditions.tsv
similarity index 100%
rename from tests/cli/expected_output/model/conditions.tsv
rename to test/cli/expected_output/model/conditions.tsv
diff --git a/tests/cli/expected_output/model/measurements.tsv b/test/cli/expected_output/model/measurements.tsv
similarity index 100%
rename from tests/cli/expected_output/model/measurements.tsv
rename to test/cli/expected_output/model/measurements.tsv
diff --git a/tests/cli/expected_output/model/model.xml b/test/cli/expected_output/model/model.xml
similarity index 100%
rename from tests/cli/expected_output/model/model.xml
rename to test/cli/expected_output/model/model.xml
diff --git a/tests/cli/expected_output/models/model_1/observables.tsv b/test/cli/expected_output/model/observables.tsv
similarity index 53%
rename from tests/cli/expected_output/models/model_1/observables.tsv
rename to test/cli/expected_output/model/observables.tsv
index ab25269a..b244d788 100644
--- a/tests/cli/expected_output/models/model_1/observables.tsv
+++ b/test/cli/expected_output/model/observables.tsv
@@ -1,2 +1,2 @@
 observableId	observableFormula	noiseFormula
-obs_x2	x2+k4	noiseParameter1_obs_x2
+obs_x2	x2+k4+k5	noiseParameter1_obs_x2
diff --git a/tests/cli/expected_output/model/parameters.tsv b/test/cli/expected_output/model/parameters.tsv
similarity index 89%
rename from tests/cli/expected_output/model/parameters.tsv
rename to test/cli/expected_output/model/parameters.tsv
index 77f65aa7..24f7485d 100644
--- a/tests/cli/expected_output/model/parameters.tsv
+++ b/test/cli/expected_output/model/parameters.tsv
@@ -3,4 +3,5 @@ k1	k_{1}	lin	0.0	1000.0	0.2	0
 k2	k_{2}	lin	0.0	1000.0	0.1	1
 k3	k_{3}	lin	0.0	1000.0	0.0	1
 k4	k_{4}	lin	0.0	1000.0	0.0	0
+k5	k_{5}	lin	0.0	1000.0	0.0	0
 sigma_x2	\sigma_{x2}	lin	1e-05	1000.0	0.15	1
diff --git a/tests/cli/expected_output/model/problem.yaml b/test/cli/expected_output/model/problem.yaml
similarity index 100%
rename from tests/cli/expected_output/model/problem.yaml
rename to test/cli/expected_output/model/problem.yaml
diff --git a/tests/cli/expected_output/models/model_1/conditions.tsv b/test/cli/expected_output/models/model_1/conditions.tsv
similarity index 100%
rename from tests/cli/expected_output/models/model_1/conditions.tsv
rename to test/cli/expected_output/models/model_1/conditions.tsv
diff --git a/tests/cli/expected_output/models/model_1/measurements.tsv b/test/cli/expected_output/models/model_1/measurements.tsv
similarity index 100%
rename from tests/cli/expected_output/models/model_1/measurements.tsv
rename to test/cli/expected_output/models/model_1/measurements.tsv
diff --git a/tests/cli/expected_output/models/model_1/model.xml b/test/cli/expected_output/models/model_1/model.xml
similarity index 100%
rename from tests/cli/expected_output/models/model_1/model.xml
rename to test/cli/expected_output/models/model_1/model.xml
diff --git a/tests/cli/expected_output/models/model_2/observables.tsv b/test/cli/expected_output/models/model_1/observables.tsv
similarity index 53%
rename from tests/cli/expected_output/models/model_2/observables.tsv
rename to test/cli/expected_output/models/model_1/observables.tsv
index ab25269a..b244d788 100644
--- a/tests/cli/expected_output/models/model_2/observables.tsv
+++ b/test/cli/expected_output/models/model_1/observables.tsv
@@ -1,2 +1,2 @@
 observableId	observableFormula	noiseFormula
-obs_x2	x2+k4	noiseParameter1_obs_x2
+obs_x2	x2+k4+k5	noiseParameter1_obs_x2
diff --git a/tests/cli/expected_output/models/model_1/parameters.tsv b/test/cli/expected_output/models/model_1/parameters.tsv
similarity index 89%
rename from tests/cli/expected_output/models/model_1/parameters.tsv
rename to test/cli/expected_output/models/model_1/parameters.tsv
index 77f65aa7..24f7485d 100644
--- a/tests/cli/expected_output/models/model_1/parameters.tsv
+++ b/test/cli/expected_output/models/model_1/parameters.tsv
@@ -3,4 +3,5 @@ k1	k_{1}	lin	0.0	1000.0	0.2	0
 k2	k_{2}	lin	0.0	1000.0	0.1	1
 k3	k_{3}	lin	0.0	1000.0	0.0	1
 k4	k_{4}	lin	0.0	1000.0	0.0	0
+k5	k_{5}	lin	0.0	1000.0	0.0	0
 sigma_x2	\sigma_{x2}	lin	1e-05	1000.0	0.15	1
diff --git a/tests/cli/expected_output/models/model_1/problem.yaml b/test/cli/expected_output/models/model_1/problem.yaml
similarity index 100%
rename from tests/cli/expected_output/models/model_1/problem.yaml
rename to test/cli/expected_output/models/model_1/problem.yaml
diff --git a/tests/cli/expected_output/models/model_2/conditions.tsv b/test/cli/expected_output/models/model_2/conditions.tsv
similarity index 100%
rename from tests/cli/expected_output/models/model_2/conditions.tsv
rename to test/cli/expected_output/models/model_2/conditions.tsv
diff --git a/tests/cli/expected_output/models/model_2/measurements.tsv b/test/cli/expected_output/models/model_2/measurements.tsv
similarity index 100%
rename from tests/cli/expected_output/models/model_2/measurements.tsv
rename to test/cli/expected_output/models/model_2/measurements.tsv
diff --git a/tests/cli/expected_output/models/model_2/model.xml b/test/cli/expected_output/models/model_2/model.xml
similarity index 100%
rename from tests/cli/expected_output/models/model_2/model.xml
rename to test/cli/expected_output/models/model_2/model.xml
diff --git a/tests/model/expected_output/petab/observables.tsv b/test/cli/expected_output/models/model_2/observables.tsv
similarity index 53%
rename from tests/model/expected_output/petab/observables.tsv
rename to test/cli/expected_output/models/model_2/observables.tsv
index ab25269a..b244d788 100644
--- a/tests/model/expected_output/petab/observables.tsv
+++ b/test/cli/expected_output/models/model_2/observables.tsv
@@ -1,2 +1,2 @@
 observableId	observableFormula	noiseFormula
-obs_x2	x2+k4	noiseParameter1_obs_x2
+obs_x2	x2+k4+k5	noiseParameter1_obs_x2
diff --git a/tests/cli/expected_output/models/model_2/parameters.tsv b/test/cli/expected_output/models/model_2/parameters.tsv
similarity index 89%
rename from tests/cli/expected_output/models/model_2/parameters.tsv
rename to test/cli/expected_output/models/model_2/parameters.tsv
index 48169160..6162eedb 100644
--- a/tests/cli/expected_output/models/model_2/parameters.tsv
+++ b/test/cli/expected_output/models/model_2/parameters.tsv
@@ -3,4 +3,5 @@ k1	k_{1}	lin	0.0	1000.0	0.2	1
 k2	k_{2}	lin	0.0	1000.0	0.1	1
 k3	k_{3}	lin	0.0	1000.0	0.0	0
 k4	k_{4}	lin	0.0	1000.0	0.0	0
+k5	k_{5}	lin	0.0	1000.0	0.0	0
 sigma_x2	\sigma_{x2}	lin	1e-05	1000.0	0.15	1
diff --git a/tests/cli/expected_output/models/model_2/problem.yaml b/test/cli/expected_output/models/model_2/problem.yaml
similarity index 100%
rename from tests/cli/expected_output/models/model_2/problem.yaml
rename to test/cli/expected_output/models/model_2/problem.yaml
diff --git a/tests/cli/input/model.yaml b/test/cli/input/model.yaml
similarity index 100%
rename from tests/cli/input/model.yaml
rename to test/cli/input/model.yaml
diff --git a/tests/cli/input/models.yaml b/test/cli/input/models.yaml
similarity index 100%
rename from tests/cli/input/models.yaml
rename to test/cli/input/models.yaml
diff --git a/tests/cli/test_cli.py b/test/cli/test_cli.py
similarity index 74%
rename from tests/cli/test_cli.py
rename to test/cli/test_cli.py
index 81c8fa67..c01982d3 100644
--- a/tests/cli/test_cli.py
+++ b/test/cli/test_cli.py
@@ -1,30 +1,33 @@
 import filecmp
 from pathlib import Path
 
+import pytest
 from click.testing import CliRunner
-#from petab_select import Model
+
+# from petab_select import Model
 import petab_select.cli
-import pytest
+
+base_dir = Path(__file__).parent
 
 
 @pytest.fixture
 def output_path() -> Path:
-    return Path('output')
+    return base_dir / 'output'
 
 
 @pytest.fixture
 def expected_output_path() -> Path:
-    return Path('expected_output')
+    return base_dir / 'expected_output'
 
 
 @pytest.fixture
-def model_yaml() -> str:
-    return 'input/model.yaml'
+def model_yaml() -> Path:
+    return base_dir / 'input' / 'model.yaml'
 
 
 @pytest.fixture
-def models_yaml() -> str:
-    return 'input/models.yaml'
+def models_yaml() -> Path:
+    return base_dir / 'input' / 'models.yaml'
 
 
 @pytest.fixture
@@ -42,14 +45,21 @@ def test_model_to_petab(
     output_path_model = output_path / 'model'
     output_path_model.mkdir(parents=True, exist_ok=True)
 
-    result = cli_runner.invoke(petab_select.cli.model_to_petab, [
-        '-y', model_yaml,
-        '-o', output_path_model,
-    ])
+    result = cli_runner.invoke(
+        petab_select.cli.model_to_petab,
+        [
+            '--model',
+            model_yaml,
+            '--output',
+            output_path_model,
+        ],
+    )
 
     print(result.stdout)
     # The new PEtab problem YAML file is output to stdout correctly.
-    assert result.stdout == 'output/model/problem.yaml\n'
+    assert (
+        result.stdout == f'{base_dir / "output" / "model" / "problem.yaml"}\n'
+    )
 
     comparison = filecmp.dircmp(
         expected_output_path / 'model',
@@ -77,16 +87,21 @@ def test_models_to_petab(
     output_path_models = output_path / 'models'
     output_path_models.mkdir(parents=True, exist_ok=True)
 
-    result = cli_runner.invoke(petab_select.cli.models_to_petab, [
-        '-y', models_yaml,
-        '-o', output_path_models,
-    ])
+    result = cli_runner.invoke(
+        petab_select.cli.models_to_petab,
+        [
+            '--models',
+            models_yaml,
+            '--output',
+            output_path_models,
+        ],
+    )
 
     # The new PEtab problem YAML files are output with model IDs to `stdout`
     # correctly.
     assert result.stdout == (
-        'model_1\toutput/models/model_1/problem.yaml\n'
-        'model_2\toutput/models/model_2/problem.yaml\n'
+        f'model_1\t{base_dir / "output" / "models" / "model_1" / "problem.yaml"}\n'
+        f'model_2\t{base_dir / "output" / "models" / "model_2" / "problem.yaml"}\n'
     )
 
     comparison = filecmp.dircmp(
diff --git a/test/model/__init__.py b/test/model/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/model/expected_output/petab/conditions.tsv b/test/model/expected_output/petab/conditions.tsv
similarity index 100%
rename from tests/model/expected_output/petab/conditions.tsv
rename to test/model/expected_output/petab/conditions.tsv
diff --git a/tests/model/expected_output/petab/measurements.tsv b/test/model/expected_output/petab/measurements.tsv
similarity index 100%
rename from tests/model/expected_output/petab/measurements.tsv
rename to test/model/expected_output/petab/measurements.tsv
diff --git a/tests/model/expected_output/petab/model.xml b/test/model/expected_output/petab/model.xml
similarity index 100%
rename from tests/model/expected_output/petab/model.xml
rename to test/model/expected_output/petab/model.xml
diff --git a/tests/cli/expected_output/model/observables.tsv b/test/model/expected_output/petab/observables.tsv
similarity index 53%
rename from tests/cli/expected_output/model/observables.tsv
rename to test/model/expected_output/petab/observables.tsv
index ab25269a..b244d788 100644
--- a/tests/cli/expected_output/model/observables.tsv
+++ b/test/model/expected_output/petab/observables.tsv
@@ -1,2 +1,2 @@
 observableId	observableFormula	noiseFormula
-obs_x2	x2+k4	noiseParameter1_obs_x2
+obs_x2	x2+k4+k5	noiseParameter1_obs_x2
diff --git a/tests/model/expected_output/petab/parameters.tsv b/test/model/expected_output/petab/parameters.tsv
similarity index 89%
rename from tests/model/expected_output/petab/parameters.tsv
rename to test/model/expected_output/petab/parameters.tsv
index 77f65aa7..24f7485d 100644
--- a/tests/model/expected_output/petab/parameters.tsv
+++ b/test/model/expected_output/petab/parameters.tsv
@@ -3,4 +3,5 @@ k1	k_{1}	lin	0.0	1000.0	0.2	0
 k2	k_{2}	lin	0.0	1000.0	0.1	1
 k3	k_{3}	lin	0.0	1000.0	0.0	1
 k4	k_{4}	lin	0.0	1000.0	0.0	0
+k5	k_{5}	lin	0.0	1000.0	0.0	0
 sigma_x2	\sigma_{x2}	lin	1e-05	1000.0	0.15	1
diff --git a/tests/model/expected_output/petab/problem.yaml b/test/model/expected_output/petab/problem.yaml
similarity index 100%
rename from tests/model/expected_output/petab/problem.yaml
rename to test/model/expected_output/petab/problem.yaml
diff --git a/tests/model/input/model.yaml b/test/model/input/model.yaml
similarity index 100%
rename from tests/model/input/model.yaml
rename to test/model/input/model.yaml
diff --git a/tests/model/test_model.py b/test/model/test_model.py
similarity index 78%
rename from tests/model/test_model.py
rename to test/model/test_model.py
index 55e8ebd5..1717b4d2 100644
--- a/tests/model/test_model.py
+++ b/test/model/test_model.py
@@ -1,21 +1,26 @@
 import filecmp
 from pathlib import Path
 
-from petab_select import Model
 import pytest
 
+from petab_select import Model
+
+base_dir = Path(__file__).parent
+
 
 @pytest.fixture
 def output_path() -> Path:
-    return Path('output')
+    return base_dir / 'output'
+
 
 @pytest.fixture
 def expected_output_path() -> Path:
-    return Path('expected_output')
+    return base_dir / 'expected_output'
+
 
 @pytest.fixture
 def model() -> Model:
-    return Model.from_yaml('input/model.yaml')
+    return Model.from_yaml(base_dir / 'input' / 'model.yaml')
 
 
 def test_model_to_petab(model, output_path, expected_output_path) -> None:
@@ -24,8 +29,9 @@ def test_model_to_petab(model, output_path, expected_output_path) -> None:
     output_path_petab.mkdir(parents=True, exist_ok=True)
     # TODO test `petab_problem`? Shouldn't be necessary since the generated
     # files are tested below.
-    petab_problem, petab_problem_yaml = \
-        model.to_petab(output_path=output_path_petab)
+    petab_problem, petab_problem_yaml = model.to_petab(
+        output_path=output_path_petab
+    )
 
     comparison = filecmp.dircmp(
         expected_output_path / 'petab',
diff --git a/test/model_space/__init__.py b/test/model_space/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/model_space/model_space_file_1.tsv b/test/model_space/model_space_file_1.tsv
similarity index 100%
rename from tests/model_space/model_space_file_1.tsv
rename to test/model_space/model_space_file_1.tsv
diff --git a/tests/model_space/model_space_file_2.tsv b/test/model_space/model_space_file_2.tsv
similarity index 100%
rename from tests/model_space/model_space_file_2.tsv
rename to test/model_space/model_space_file_2.tsv
diff --git a/tests/model_space/test_model_space.py b/test/model_space/test_model_space.py
similarity index 65%
rename from tests/model_space/test_model_space.py
rename to test/model_space/test_model_space.py
index 2b0e2f99..de690713 100644
--- a/tests/model_space/test_model_space.py
+++ b/test/model_space/test_model_space.py
@@ -1,35 +1,32 @@
+from pathlib import Path
 from typing import List
 
 import pandas as pd
 import pytest
 
-from petab_select.constants import (
-    ESTIMATE,
-
-    MODEL_SUBSPACE_ID,
-    PETAB_YAML,
-    PARAMETER_VALUE_DELIMITER,
-)
-from petab_select.model import (
-    Model,
-)
-from petab_select.model_space import (
-    ModelSpace,
-)
-
 from petab_select.candidate_space import (
     BackwardCandidateSpace,
     BruteForceCandidateSpace,
     ForwardCandidateSpace,
     LateralCandidateSpace,
 )
+from petab_select.constants import (
+    ESTIMATE,
+    MODEL_SUBSPACE_ID,
+    PARAMETER_VALUE_DELIMITER,
+    PETAB_YAML,
+)
+from petab_select.model import Model
+from petab_select.model_space import ModelSpace
+
+base_dir = Path(__file__).parent
 
 
 @pytest.fixture
-def model_space_files() -> List[str]:
+def model_space_files() -> List[Path]:
     return [
-        'model_space_file_1.tsv',
-        'model_space_file_2.tsv',
+        base_dir / 'model_space_file_1.tsv',
+        base_dir / 'model_space_file_2.tsv',
     ]
 
 
@@ -46,9 +43,18 @@ def test_model_space_forward_virtual(model_space):
     # model is used. This means the expected models are the "smallest" models (as many
     # fixed parameters as possible) in the model space.
     expected_models = [
-        ('model_subspace_1', {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.0}),
-        ('model_subspace_1', {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.1}),
-        ('model_subspace_2', {'k1': 0.0, 'k2': 0.0, 'k3': 0.0, 'k4': ESTIMATE}),
+        (
+            'model_subspace_1',
+            {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.0},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.1},
+        ),
+        (
+            'model_subspace_2',
+            {'k1': 0.0, 'k2': 0.0, 'k3': 0.0, 'k4': ESTIMATE},
+        ),
     ]
 
     models = [
@@ -57,15 +63,9 @@ def test_model_space_forward_virtual(model_space):
     ]
 
     # Search found only expected models.
-    assert all([
-        model in expected_models
-        for model in models
-    ])
+    assert all([model in expected_models for model in models])
     # All expected models have now been added to the candidate space.
-    assert all([
-        model in models
-        for model in expected_models
-    ])
+    assert all([model in models for model in expected_models])
     # Probably unnecessary: same number of models in expectation vs realization
     assert len(expected_models) == len(candidate_space.models)
 
@@ -79,11 +79,11 @@ def test_model_space_backward_virtual(model_space):
     # model is used. This means the expected models are the "smallest" models (as many
     # fixed parameters as possible) in the model space.
     expected_models = [
-        ('model_subspace_1', {f'k{i}': ESTIMATE for i in range(1,5)}),
+        ('model_subspace_1', {f'k{i}': ESTIMATE for i in range(1, 5)}),
         # This model is not included because it is exactly the same as the
         # other model (same PEtab YAML and parameterization), hence has been
         # excluded from the candidate space.
-        #('model_subspace_3', {f'k{i}': ESTIMATE for i in range(1,5)}),
+        # ('model_subspace_3', {f'k{i}': ESTIMATE for i in range(1,5)}),
     ]
 
     models = [
@@ -92,15 +92,9 @@ def test_model_space_backward_virtual(model_space):
     ]
 
     # Search found only expected models.
-    assert all([
-        model in expected_models
-        for model in models
-    ])
+    assert all([model in expected_models for model in models])
     # All expected models have now been added to the candidate space.
-    assert all([
-        model in models
-        for model in expected_models
-    ])
+    assert all([model in models for model in expected_models])
     # Probably unnecessary: same number of models in expectation vs realization
     assert len(expected_models) == len(candidate_space.models)
 
@@ -112,19 +106,58 @@ def test_model_space_brute_force_limit(model_space):
     # There are fifteen total models in the model space. Limiting to 13 models should
     # result in all models except the last two models in the last model subspace.
     expected_models = [
-        ('model_subspace_1', {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.0}),
-        ('model_subspace_1', {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.1}),
-        ('model_subspace_1', {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': ESTIMATE}),
-        ('model_subspace_1', {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.0}),
-        ('model_subspace_1', {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.1}),
-        ('model_subspace_1', {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': ESTIMATE}),
-        ('model_subspace_1', {'k1': ESTIMATE, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.0}),
-        ('model_subspace_1', {'k1': ESTIMATE, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.1}),
-        ('model_subspace_1', {'k1': ESTIMATE, 'k2': 0.1, 'k3': ESTIMATE, 'k4': ESTIMATE}),
-        ('model_subspace_1', {'k1': ESTIMATE, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.0}),
-        ('model_subspace_1', {'k1': ESTIMATE, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.1}),
-        ('model_subspace_1', {'k1': ESTIMATE, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': ESTIMATE}),
-        ('model_subspace_2', {'k1': 0.0, 'k2': 0.0, 'k3': 0.0, 'k4': ESTIMATE}),
+        (
+            'model_subspace_1',
+            {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.0},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.1},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': ESTIMATE},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.0},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.1},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': ESTIMATE},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': ESTIMATE, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.0},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': ESTIMATE, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.1},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': ESTIMATE, 'k2': 0.1, 'k3': ESTIMATE, 'k4': ESTIMATE},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': ESTIMATE, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.0},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': ESTIMATE, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.1},
+        ),
+        (
+            'model_subspace_1',
+            {'k1': ESTIMATE, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': ESTIMATE},
+        ),
+        (
+            'model_subspace_2',
+            {'k1': 0.0, 'k2': 0.0, 'k3': 0.0, 'k4': ESTIMATE},
+        ),
     ]
 
     models = [
@@ -133,18 +166,13 @@ def test_model_space_brute_force_limit(model_space):
     ]
 
     # Search found only expected models.
-    assert all([
-        model in expected_models
-        for model in models
-    ])
+    assert all([model in expected_models for model in models])
     # All expected models have now been added to the candidate space.
-    assert all([
-        model in models
-        for model in expected_models
-    ])
+    assert all([model in models for model in expected_models])
     # Probably unnecessary: same number of models in expectation vs realization
     assert len(expected_models) == len(candidate_space.models)
 
+
 '''
 @pytest.fixture
 def e():
diff --git a/test/model_subspace/__init__.py b/test/model_subspace/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/model_subspace/test_model_subspace.py b/test/model_subspace/test_model_subspace.py
similarity index 56%
rename from tests/model_subspace/test_model_subspace.py
rename to test/model_subspace/test_model_subspace.py
index 60e4815f..71dbff04 100644
--- a/tests/model_subspace/test_model_subspace.py
+++ b/test/model_subspace/test_model_subspace.py
@@ -1,14 +1,15 @@
 from pathlib import Path
 
-from more_itertools import one
 import pandas as pd
 import pytest
+from more_itertools import one
 
 import petab_select
 from petab_select.candidate_space import (
     BackwardCandidateSpace,
     BruteForceCandidateSpace,
     ForwardCandidateSpace,
+    LateralCandidateSpace,
 )
 from petab_select.constants import (
     ESTIMATE,
@@ -16,38 +17,49 @@
     PARAMETER_VALUE_DELIMITER,
     PETAB_YAML,
 )
-from petab_select.model import (
-    Model,
-)
-from petab_select.model_subspace import (
-    ModelSubspace,
-)
+from petab_select.model import Model
+from petab_select.model_subspace import ModelSubspace
 
 
 @pytest.fixture
-def definition() -> pd.Series:
+def model_subspace_id_and_definition() -> pd.Series:
+    model_subspace_id = 'model_subspace_1'
     data = {
-        MODEL_SUBSPACE_ID: 'model_subspace_1',
-        PETAB_YAML: Path(__file__).parent.parent.parent / 'doc' / 'examples' / 'model_selection' / 'petab_problem.yaml',  # noqa: E501
+        PETAB_YAML: Path(__file__).parent.parent.parent
+        / 'doc'
+        / 'examples'
+        / 'model_selection'
+        / 'petab_problem.yaml',
         'k1': 0.2,
         'k2': PARAMETER_VALUE_DELIMITER.join(['0.1', ESTIMATE]),
         'k3': ESTIMATE,
         'k4': PARAMETER_VALUE_DELIMITER.join(['0', '0.1', ESTIMATE]),
     }
-    return pd.Series(data=data, dtype=str)
+    return model_subspace_id, pd.Series(data=data, dtype=str)
 
 
 @pytest.fixture
-def model_subspace(definition) -> ModelSubspace:
-    return petab_select.model_subspace.ModelSubspace.from_definition(definition)
+def model_subspace(model_subspace_id_and_definition) -> ModelSubspace:
+    model_subspace_id, definition = model_subspace_id_and_definition
+    return petab_select.model_subspace.ModelSubspace.from_definition(
+        model_subspace_id=model_subspace_id,
+        definition=definition,
+    )
 
 
 @pytest.fixture
 def initial_model(model_subspace) -> Model:
     estimated_parameters = ['k3', 'k4']
-    model = one(model_subspace.get_models(estimated_parameters=estimated_parameters))
+    model = one(
+        model_subspace.get_models(estimated_parameters=estimated_parameters)
+    )
     # Initial model is parameterized as expected.
-    assert model.parameters == {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': ESTIMATE}
+    assert model.parameters == {
+        'k1': 0.2,
+        'k2': 0.1,
+        'k3': ESTIMATE,
+        'k4': ESTIMATE,
+    }
     return model
 
 
@@ -57,7 +69,11 @@ def test_from_definition(model_subspace):
     assert model_subspace.model_subspace_id == 'model_subspace_1'
     # PEtab YAML is parsed
     assert model_subspace.petab_yaml.samefile(
-        Path(__file__).parent.parent.parent / 'doc' / 'examples' / 'model_selection' / 'petab_problem.yaml',  # noqa: E501
+        Path(__file__).parent.parent.parent
+        / 'doc'
+        / 'examples'
+        / 'model_selection'
+        / 'petab_problem.yaml',
     )
     # Fixed parameters are parsed
     assert model_subspace.parameters['k1'] == [0.2]
@@ -73,22 +89,28 @@ def test_from_definition(model_subspace):
 def test_get_models(model_subspace):
     """The getter for models with specific estimated parameters works."""
     estimated_parameters = ['k2', 'k3']
-    models = list(model_subspace.get_models(estimated_parameters=estimated_parameters))
-    expected_model_parameterizations = [
+    models = list(
+        model_subspace.get_models(estimated_parameters=estimated_parameters)
+    )
+    expected_parameterizations = [
         {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.0},
         {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.1},
     ]
-    model_parameterizations = [model.parameters for model in models]
+    test_parameterizations = [model.parameters for model in models]
     # Getter gets only expected models.
-    assert all([
-        parameterization in expected_model_parameterizations
-        for parameterization in model_parameterizations
-    ])
+    assert all(
+        [
+            test_parameterization in expected_parameterizations
+            for test_parameterization in expected_parameterizations
+        ]
+    )
     # Getter gets all expected models.
-    assert all([
-        parameterization in model_parameterizations
-        for parameterization in expected_model_parameterizations
-    ])
+    assert all(
+        [
+            expected_parameterization in test_parameterizations
+            for expected_parameterization in expected_parameterizations
+        ]
+    )
 
 
 def test_search_forward(model_subspace, initial_model):
@@ -99,9 +121,13 @@ def test_search_forward(model_subspace, initial_model):
     # Only one model is possible in the forward direction.
     assert len(candidate_space.models) == 1
     # The one model has the expected parameterization.
-    expected_forward_parameterization = \
-        {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': ESTIMATE}
-    assert one(candidate_space.models).parameters == expected_forward_parameterization
+    expected_parameterization = {
+        'k1': 0.2,
+        'k2': ESTIMATE,
+        'k3': ESTIMATE,
+        'k4': ESTIMATE,
+    }
+    assert one(candidate_space.models).parameters == expected_parameterization
 
     # Test limit via model subspace
     # FIXME currently only returns 1 model anyway
@@ -112,8 +138,8 @@ def test_search_forward(model_subspace, initial_model):
         candidate_space=candidate_space,
         limit=limit_considered_candidates,
     )
-    # Test limit: the number of candidate models is at the limit (all models in this
-    # case).
+    # Test limit: the number of candidate models is at the limit (all models in
+    # this case).
     assert len(candidate_space.models) == limit_considered_candidates
 
 
@@ -124,24 +150,27 @@ def test_search_backward(model_subspace, initial_model):
     model_subspace.search(candidate_space=candidate_space)
     # Only two models are possible in the backward direction.
     assert len(candidate_space.models) == 2
-    expected_backward_parameterizations = [
+    expected_parameterizations = [
         {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0},
         {'k1': 0.2, 'k2': 0.1, 'k3': ESTIMATE, 'k4': 0.1},
     ]
-    model_parameterizations = [
-        model.parameters
-        for model in candidate_space.models
+    test_parameterizations = [
+        model.parameters for model in candidate_space.models
     ]
     # Search found only expected models.
-    assert all([
-        parameterization in expected_backward_parameterizations
-        for parameterization in model_parameterizations
-    ])
+    assert all(
+        [
+            test_parameterization in expected_parameterizations
+            for test_parameterization in test_parameterizations
+        ]
+    )
     # Search found all expected models.
-    assert all([
-        parameterization in model_parameterizations
-        for parameterization in expected_backward_parameterizations
-    ])
+    assert all(
+        [
+            expected_parameterization in test_parameterizations
+            for expected_parameterization in expected_parameterizations
+        ]
+    )
 
     # Test limit via model subspace
     limit_considered_candidates = 1
@@ -172,20 +201,23 @@ def test_search_brute_force(model_subspace):
         {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.1},
         {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': ESTIMATE},
     ]
-    parameterizations = [
-        model.parameters
-        for model in candidate_space.models
+    test_parameterizations = [
+        model.parameters for model in candidate_space.models
     ]
     # Search found only expected models.
-    assert all([
-        parameterization in expected_parameterizations
-        for parameterization in parameterizations
-    ])
+    assert all(
+        [
+            test_parameterization in expected_parameterizations
+            for test_parameterization in test_parameterizations
+        ]
+    )
     # Search found all expected models.
-    assert all([
-        parameterization in parameterizations
-        for parameterization in parameterizations
-    ])
+    assert all(
+        [
+            expected_parameterization in test_parameterizations
+            for expected_parameterization in expected_parameterizations
+        ]
+    )
 
     limit_accepted_candidates = 3
     candidate_space.reset(limit=limit_accepted_candidates)
@@ -199,36 +231,41 @@ def test_search_brute_force(model_subspace):
     """
     # Test limit: the number of candidate models is at the limit.
     assert len(candidate_space.models) == limit_accepted_candidates
-    # Test limit: a warning is emitted, as the candidate space is already at its limit
-    # of models.
-    with pytest.warns(RuntimeWarning, match="The candidate space has already reached its limit of accepted models.") as warning_record:
-        model_subspace.search(candidate_space=candidate_space)
-    # The search stopped after the limit was reached, hence only warned once.
-    assert len(warning_record) == 1
+    model_subspace.search(candidate_space=candidate_space)
     # Test limit: the search adds no additional models to the "full" candidate space.
     assert len(candidate_space.models) == limit_accepted_candidates
 
     limit_accepted_candidates = 6
     candidate_space.limit.set_limit(limit_accepted_candidates)
-    with pytest.warns(RuntimeWarning, match="Model has been previously excluded from the candidate space so is skipped here.") as warning_record:
+    with pytest.warns(
+        RuntimeWarning,
+        match="Model has been previously excluded from the candidate space so is skipped here.",
+    ) as warning_record:
         model_subspace.search(candidate_space=candidate_space)
+    test_parameterizations = [
+        model.parameters for model in candidate_space.models
+    ]
     # Three models were excluded from the candidate space in the previous code block.
     assert len(warning_record) == 3
     # Test limit: the number of candidate models is at the limit (all models in this
     # case).
     assert len(candidate_space.models) == limit_accepted_candidates
     # Search found only expected models.
-    assert all([
-        parameterization in expected_parameterizations
-        for parameterization in parameterizations
-    ])
+    assert all(
+        [
+            test_parameterization in expected_parameterizations
+            for test_parameterization in test_parameterizations
+        ]
+    )
     # Test exclusions: all models have now been added to the candidate space.
     # TODO ideally with only 3 additional calls to `candidate_space.consider`, assuming
     #      the model_subspace excluded the first three models it had already sent.
-    assert all([
-        parameterization in parameterizations
-        for parameterization in expected_parameterizations
-    ])
+    assert all(
+        [
+            expected_parameterization in test_parameterizations
+            for expected_parameterization in expected_parameterizations
+        ]
+    )
 
     # Test limit via model subspace
     limit_considered_candidates = 1
@@ -241,3 +278,47 @@ def test_search_brute_force(model_subspace):
     # Test limit: the number of candidate models is at the limit (all models in this
     # case).
     assert len(candidate_space.models) == limit_considered_candidates
+
+
+def test_search_swap(model_subspace, initial_model):
+    # TODO exclude history, use limit
+    candidate_space = LateralCandidateSpace(predecessor_model=initial_model)
+
+    model_subspace.search(candidate_space=candidate_space)
+    # Only two models are possible in the swap direction.
+    assert len(candidate_space.models) == 2
+    # The two models have the expected parameterization.
+    expected_parameterizations = [
+        {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0},
+        {'k1': 0.2, 'k2': ESTIMATE, 'k3': ESTIMATE, 'k4': 0.1},
+    ]
+    test_parameterizations = [
+        model.parameters for model in candidate_space.models
+    ]
+    # Search found only expected models.
+    assert all(
+        [
+            test_parameterization in expected_parameterizations
+            for test_parameterization in test_parameterizations
+        ]
+    )
+    # Search found all expected models.
+    assert all(
+        [
+            expected_parameterization in test_parameterizations
+            for expected_parameterization in expected_parameterizations
+        ]
+    )
+
+    # Test limit via model subspace
+    # FIXME currently only returns 1 model anyway
+    limit_considered_candidates = 1
+    model_subspace.reset_exclusions()
+    candidate_space.reset(predecessor_model=initial_model)
+    model_subspace.search(
+        candidate_space=candidate_space,
+        limit=limit_considered_candidates,
+    )
+    # Test limit: the number of candidate models is at the limit (all models in this
+    # case).
+    assert len(candidate_space.models) == limit_considered_candidates
diff --git a/test/pypesto/__init__.py b/test/pypesto/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/pypesto/generate_expected_models.py b/test/pypesto/generate_expected_models.py
similarity index 73%
rename from tests/pypesto/generate_expected_models.py
rename to test/pypesto/generate_expected_models.py
index d1bab3b7..2d04d815 100644
--- a/tests/pypesto/generate_expected_models.py
+++ b/test/pypesto/generate_expected_models.py
@@ -1,19 +1,15 @@
 from pathlib import Path
 
 import fides
-from more_itertools import one
 import pandas as pd
-import petab_select
-from petab_select import Model
-from petab_select.constants import (
-    CRITERIA,
-    ESTIMATED_PARAMETERS,
-    MODEL,
-)
 import pypesto.engine
 import pypesto.optimize
 import pypesto.select
+from more_itertools import one
 
+import petab_select
+from petab_select import Model
+from petab_select.constants import CRITERIA, ESTIMATED_PARAMETERS, MODEL
 
 SKIP_TEST_CASES_WITH_PREEXISTING_EXPECTED_MODEL = False
 
@@ -28,8 +24,9 @@
 # Reduce runtime but with high reproducibility
 minimize_options = {
     'n_starts': 100,
-    'optimizer': pypesto.optimize.FidesOptimizer(),
+    'optimizer': pypesto.optimize.FidesOptimizer(verbose=0),
     'engine': pypesto.engine.MultiProcessEngine(),
+    'filename': None,
 }
 
 # Indentation to match `test_pypesto.py`, to make it easier to keep files similar.
@@ -52,17 +49,21 @@
         petab_select_problem = petab_select.Problem.from_yaml(
             test_case_path / 'petab_select_problem.yaml',
         )
-        pypesto_select_problem = \
-            pypesto.select.Problem(petab_select_problem=petab_select_problem)
+        pypesto_select_problem = pypesto.select.Problem(
+            petab_select_problem=petab_select_problem
+        )
 
         # Run the selection process until "exhausted".
         pypesto_select_problem.select_to_completion(
             minimize_options=minimize_options,
         )
 
-        # Get the best model, load the expected model.
-        models = pypesto_select_problem.history.values()
-        best_model = petab_select_problem.get_best(models)
+        # Get the best model
+        best_model = petab_select_problem.get_best(
+            models=pypesto_select_problem.calibrated_models.values(),
+        )
 
         # Generate the expected model.
-        best_model.to_yaml(expected_model_yaml, paths_relative_to=test_case_path)
+        best_model.to_yaml(
+            expected_model_yaml, paths_relative_to=test_case_path
+        )
diff --git a/tests/pypesto/test_pypesto.py b/test/pypesto/test_pypesto.py
similarity index 80%
rename from tests/pypesto/test_pypesto.py
rename to test/pypesto/test_pypesto.py
index 9234cbfc..5c333002 100644
--- a/tests/pypesto/test_pypesto.py
+++ b/test/pypesto/test_pypesto.py
@@ -1,21 +1,17 @@
 from pathlib import Path
 
 import fides
-from more_itertools import one
 import numpy as np
 import pandas as pd
-import petab_select
-from petab_select import Model
-from petab_select.constants import (
-    CRITERIA,
-    ESTIMATED_PARAMETERS,
-    MODEL,
-)
 import pypesto.engine
 import pypesto.optimize
 import pypesto.select
 import pytest
+from more_itertools import one
 
+import petab_select
+from petab_select import Model
+from petab_select.constants import CRITERIA, ESTIMATED_PARAMETERS, MODEL
 
 # Set to `[]` to test all
 test_cases = [
@@ -28,10 +24,12 @@
 # Reduce runtime but with high reproducibility
 minimize_options = {
     'n_starts': 10,
-    'optimizer': pypesto.optimize.FidesOptimizer(),
+    'optimizer': pypesto.optimize.FidesOptimizer(verbose=0),
     'engine': pypesto.engine.MultiProcessEngine(),
+    'filename': None,
 }
 
+
 def test_pypesto():
     for test_case_path in test_cases_path.glob('*'):
         if test_cases and test_case_path.stem not in test_cases:
@@ -42,8 +40,9 @@ def test_pypesto():
         petab_select_problem = petab_select.Problem.from_yaml(
             test_case_path / 'petab_select_problem.yaml',
         )
-        pypesto_select_problem = \
-            pypesto.select.Problem(petab_select_problem=petab_select_problem)
+        pypesto_select_problem = pypesto.select.Problem(
+            petab_select_problem=petab_select_problem
+        )
 
         # Run the selection process until "exhausted".
         pypesto_select_problem.select_to_completion(
@@ -51,8 +50,9 @@ def test_pypesto():
         )
 
         # Get the best model, load the expected model.
-        models = pypesto_select_problem.history.values()
-        best_model = petab_select_problem.get_best(models)
+        best_model = petab_select_problem.get_best(
+            models=pypesto_select_problem.calibrated_models.values(),
+        )
         expected_model = Model.from_yaml(expected_model_yaml)
 
         def get_series(model) -> pd.Series:
@@ -60,6 +60,7 @@ def get_series(model) -> pd.Series:
                 getattr(model, dict_attribute),
                 dtype=np.float64,
             ).sort_index()
+
         # The estimated parameters and criteria values are as expected.
         for dict_attribute in [CRITERIA, ESTIMATED_PARAMETERS]:
             pd.testing.assert_series_equal(