Repository: popsim-consortium/stdpopsim (language: python)

Issue:
Allele frequency conditioning with multiple `DrawMutation` events at the same locus
Relates to #608.

If multiple `DrawMutation` events with the same single site id are supplied to the SLiM engine, then multiple mutations of the same type will be generated. If these mutations make it into the sampled individuals, they will ultimately end up as different alleles (integer codes) in the outputted tree sequence.

This creates a bit of a dilemma for `ConditionOnAlleleFrequency`. What exactly should be conditioned on? Some options are:

1. The simulation could abort if there's more than one `DrawMutation` event at the same site. This is *de facto* what happens currently -- if there's more than one extant mutation at a site when the AF condition is evaluated, then SLiM will exit with a rather ambiguous error message. So we'd just need to add a check for overlapping `DrawMutation` events in `slim_makescript` and exit with an informative error.
2. The allele frequency condition could apply to the sum of frequencies of derived alleles.
3. The allele frequency condition could apply independently to all derived alleles. This seems useful in some cases, like for conditioning on the survival of all derived alleles. However we'd need some additional checks to ensure that the condition is possible to satisfy.
4. Separate conditions could be applied to separate alleles. I think the only reasonable way to allow this would be to permit multiple mutation types at the same "single site" (in which case alleles could have separate fitness coefficients) -- but this would require more infrastructure, testing, etc.

My vote is for (1) as this will require very little work to implement, and should cover many use cases (?). Ultimately it'd be nice to extend this to (4).
fix allele frequency calculation in slim code, or disallow multiple drawn mutations per MutationType
The allele frequency calculation in the slim engine (used for allele frequency conditioning) currently assumes there is only one mutation per `MutationType`.
https://github.com/popsim-consortium/stdpopsim/blob/973d5b05f853f5412fceff5b48e12e3c5573eb2d/stdpopsim/slim_engine.py#L444-L452

But a user could easily break this by drawing two mutations of the same mutation type. See note here: https://github.com/popsim-consortium/stdpopsim/blob/973d5b05f853f5412fceff5b48e12e3c5573eb2d/stdpopsim/ext/selection.py#L115-L133
update ubuntu in tests
Github warns us about our tests now:
```
The ubuntu-18.04 environment is deprecated, consider switching to ubuntu-20.04(ubuntu-latest), or ubuntu-22.04 instead. For more details see https://github.com/actions/virtual-environments/issues/6002
```
... I figure we should just change to [20.04 here](https://github.com/popsim-consortium/stdpopsim/blob/main/.github/workflows/tests.yml#L20).


Discussion / hints:
I vote for (1) also!



Fix patch:
diff --git a/docs/selection_example.py b/docs/selection_example.py
index 16590b39b..9e6e0311a 100644
--- a/docs/selection_example.py
+++ b/docs/selection_example.py
@@ -4,6 +4,10 @@
 
 logger = logging.getLogger(__name__)
 
+# TODO: This example has been updated to reflect changes in the extended events
+# API (see PRs #1306 and #1341) but it should be run and checked for
+# correctness at some point
+
 
 def adaptive_introgression(seed):
     """
@@ -21,13 +25,6 @@ def adaptive_introgression(seed):
         100, 0, 0, 100, 2, 2  # YRI, CEU, CHB, Papuan, DenA, NeaA
     )
 
-    # One mutation type, which we'll use for the positively selected mutation.
-    # Neutral mutations will be added by the SLiM engine as usual, after the
-    # SLiM phase of the simulation has completed.
-    positive = stdpopsim.MutationType(convert_to_substitution=False)
-    mutation_types = [positive]
-    mut_id = len(mutation_types) - 1
-
     # We need some demographic model parameters to set bounds on the timing
     # of random variables and extended_events (below).
     # These values were copied from the PapuansOutOfAfrica_10J19 model
@@ -55,9 +52,9 @@ def adaptive_introgression(seed):
     logger.info(f"Parameters: T_mut={T_mut:.3f}, T_sel={T_sel:.3f}, s={s:.3g}")
 
     # Place the drawn mutation in the middle of the contig.
+    locus_id = "introgressed_locus"
     coordinate = round(contig.recombination_map.sequence_length / 2)
-
-    pop = {p.name: i for i, p in enumerate(model.populations)}
+    contig.add_single_site(id=locus_id, coordinate=coordinate)
 
     # Thinking forwards in time, we define a number of extended events that
     # correspond to drawing the mutation, conditioning on the new allele not
@@ -69,9 +66,8 @@ def adaptive_introgression(seed):
         # Draw mutation in DenA.
         stdpopsim.ext.DrawMutation(
             time=T_mut,
-            mutation_type_id=mut_id,
-            population_id=pop["DenA"],
-            coordinate=coordinate,
+            single_site_id=locus_id,
+            population="DenA",
         ),
         # Because the drawn mutation is neutral at the time of introduction,
         # it's likely to be lost due to drift. To avoid this, we condition on
@@ -86,8 +82,8 @@ def adaptive_introgression(seed):
             #       which will give an error due to "start_time < end_time".
             start_time=stdpopsim.ext.GenerationAfter(T_mut),
             end_time=T_Den_split,
-            mutation_type_id=mut_id,
-            population_id=pop["DenA"],
+            single_site_id=locus_id,
+            population="DenA",
             op=">",
             allele_frequency=0,
         ),
@@ -96,8 +92,8 @@ def adaptive_introgression(seed):
         stdpopsim.ext.ConditionOnAlleleFrequency(
             start_time=stdpopsim.ext.GenerationAfter(T_Den_split),
             end_time=T_mig,
-            mutation_type_id=mut_id,
-            population_id=pop["Den1"],
+            single_site_id=locus_id,
+            population="Den1",
             op=">",
             allele_frequency=0,
         ),
@@ -106,8 +102,8 @@ def adaptive_introgression(seed):
         stdpopsim.ext.ConditionOnAlleleFrequency(
             start_time=stdpopsim.ext.GenerationAfter(T_mig),
             end_time=0,
-            mutation_type_id=mut_id,
-            population_id=pop["Papuan"],
+            single_site_id=locus_id,
+            population="Papuan",
             op=">",
             allele_frequency=0,
         ),
@@ -117,8 +113,8 @@ def adaptive_introgression(seed):
         stdpopsim.ext.ChangeMutationFitness(
             start_time=T_sel,
             end_time=0,
-            mutation_type_id=mut_id,
-            population_id=pop["Papuan"],
+            single_site_id=locus_id,
+            population="Papuan",
             selection_coeff=s,
             dominance_coeff=0.5,
         ),
@@ -126,8 +122,8 @@ def adaptive_introgression(seed):
         stdpopsim.ext.ConditionOnAlleleFrequency(
             start_time=0,
             end_time=0,
-            mutation_type_id=mut_id,
-            population_id=pop["Papuan"],
+            single_site_id=locus_id,
+            population="Papuan",
             op=">",
             allele_frequency=0.05,
         ),
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
index f5c0d782f..5dcec64f7 100644
--- a/docs/tutorial.rst
+++ b/docs/tutorial.rst
@@ -1189,7 +1189,7 @@ but nonetheless there is lower diversity in exons than outside of them:
     is preliminary, and subject to change!
 
 You may be interested in simulating and tracking a single mutation. To illustrate
-this scenario, let's simulate a selective sweep until it reaches an abitrary
+this scenario, let's simulate a selective sweep until it reaches an arbitrary
 allele frequency.
 
 First, let's define a contig and a demographic model; here, we are simulating a
@@ -1212,34 +1212,21 @@ We must also decide the time the mutation will be added, when selection will
 start and at what frequency we want our selected mutation to be at the end of
 the simulation.
 
-Let's assume the mutation appeared 1000 generations ago, it has a positive
-effect on fitness (s=0.5). Also, we want the mutation to have reached a frequency
-of at least 0.8 by the end. Next, we'll walk through the steps required to do this:
-
-.. note::
-
-    Note that because we are doing a forward-in-time simulation, you should be
-    careful with your conditioning. For example, even a strongly selected mutation
-    would not be able to reach 80% frequency in just a few generations. Since
-    this conditioning works by re-running the simulation until the condition is
-    achieved, a nearly impossible condition will result in very long run times.
-
-To set things up, we need to first add the site at which the selected mutation
-will occur.  This is like adding a DFE, except to a single site -- we're saying
-that there is a potential mutation at a particular site with defined fitness
-consequences. So that we can refer to the single site later, we give it a
-unique string ID. Here, we'll add the site in the middle of the contig, with ID
-"hard sweep".
+First, we need to add the site at which the selected mutation will occur.  This
+is like adding a DFE, except to a single site -- we're saying that there is a
+potential mutation at a particular site with defined fitness consequences. So
+that we can refer to the single site later, we give it a unique string ID.
+Here, we'll add the site in the middle of the contig with ID "hard sweep",
+so named because we will imagine this beneficial mutation originates at
+frequency :math:`1 / 2N`.
 
 .. code-block:: python
 
-    mut_id = "hard sweep"
+    locus_id = "hard sweep"
     coordinate = round(contig.recombination_map.sequence_length / 2)
     contig.add_single_site(
-        id=mut_id,
+        id=locus_id,
         coordinate=coordinate,
-        fitness_coeff=0.1,
-        dominance_coeff=0.5,
     )
 
 .. note::
@@ -1251,62 +1238,34 @@ unique string ID. Here, we'll add the site in the middle of the contig, with ID
     "overwritten" and an error will be raised in simulation.
 
 Next, we will set up the "extended events" which will modify the demography.
-The first extended event is the origination of the selected mutation, which
-will occur in a random individual from the first population (id 0), 1000
-generations ago.
+This is done through :func:`stdpopsim.ext.selective_sweep`, which represents a
+general model for a mutation that is beneficial within a single population.  We
+specify that the mutation should originate 1000 generations ago in a random
+individual from the first population (named "pop_0" by default); that the
+selection coefficient for the mutation should be 0.5; and that the frequency of
+the mutation in the present day (e.g. at the end of the sweep) should be
+greater than 0.8.
 
 .. code-block:: python
 
-    T_mut = 1000
-    extended_events = [
-        stdpopsim.ext.DrawMutation(
-            time=T_mut,
-            single_site_id=mut_id,
-            population_id=0,
-        )
-    ]
-
-Next, we condition on the mutation not being lost.  Since in the next step we
-condition on the mutation being at 80% frequency at the end, this is redundant,
-but it allows the simulation to immediately restart from any generation in
-which the mutation is lost, rather than waiting until the end.  Note that this
-conditioning must start one generation after the mutation is placed, for which
-we use ``stdpopsim.ext.GenerationAfter(T_mut)``.  We cannot simply specify
-``T_mut - 1`` if rescaling is present, otherwise the conditioning would start
-at the same generation when the mutation is placed.
-
-.. code-block:: python
-
-    extended_events.append(
-        stdpopsim.ext.ConditionOnAlleleFrequency(
-            start_time=stdpopsim.ext.GenerationAfter(T_mut),
-            end_time=0,
-            single_site_id=mut_id,
-            population_id=0,
-            op=">",
-            allele_frequency=0.0,
-        )
+    extended_events = stdpopsim.ext.selective_sweep(
+        single_site_id=locus_id,
+        population="pop_0",
+        selection_coeff=0.5,
+        mutation_generation_ago=1000,
+        min_freq_at_end=0.8,
     )
 
-Finally, we condition on the mutation being above 80% at the end of the simulation.
-(The "end" is at time 0, since "time" is in generations before the end of the simulation.)
-
-.. code-block:: python
+.. note::
 
-    extended_events.append(
-        stdpopsim.ext.ConditionOnAlleleFrequency(
-            start_time=0,
-            end_time=0,
-            single_site_id=mut_id,
-            population_id=0,
-            op=">=",
-            allele_frequency=0.8,
-        )
-    )
+    Note that because we are doing a forward-in-time simulation, you should be
+    careful with your conditioning. For example, even a strongly selected mutation
+    would not be able to reach 80% frequency in just a few generations. Since
+    this conditioning works by re-running the simulation until the condition is
+    achieved, a nearly impossible condition will result in very long run times.
 
-Now we can simulate, using SLiM of course.
-For comparison, we will run the same simulation
-without selection - i.e., without the "extended events":
+Now we can simulate, using SLiM of course.  For comparison, we will run the
+same simulation without selection - i.e., without the "extended events":
 
 .. code-block:: python
 
diff --git a/stdpopsim/ext/selection.py b/stdpopsim/ext/selection.py
index 60fb6d597..46b279ff9 100644
--- a/stdpopsim/ext/selection.py
+++ b/stdpopsim/ext/selection.py
@@ -48,8 +48,9 @@ def validate_time_range(start_time, end_time):
 class ExtendedEvent(object):
     """
     ExtendedEvent is analogous to msprime.DemographicEvent, but is used here
-    for SLiM things that msprime doesn't support. Times are all in units of
-    generations before the present, just like for msprime.DemographicEvent.
+    for explicitly manipulating mutations. Currently, these are only used by
+    the SLiM engine. Times are all in units of generations before the present,
+    just like for msprime.DemographicEvent.
     """
 
     pass
@@ -58,23 +59,25 @@ class ExtendedEvent(object):
 @attr.s(kw_only=True)
 class DrawMutation(ExtendedEvent):
     """
-    TODO: docstring.
-
-    Draw a new mutation with the given mutation type in the given population
-    at the given genomic coordinate. The mutation will be added to one randomly
-    chosen chromosome in the given population.
-
-    FIXME: Drawing multiple mutations using the same mutation type causes
-           erroneous allele frequency calculation in the SLiM code!
-           Multiple drawn mutations should be disallowed, or the allele frequency
-           calculation should somehow be tied to specific drawn mutations.
-
-    See SLiM documentation for addNewDrawnMutation().
+    Introduce a new mutation with the given mutation type in the given
+    population at a given single site. The mutation will be added to one
+    randomly chosen chromosome in the given population. Only one mutation
+    may be drawn per single site.
+
+    :ivar time: The number of generations in the past at which to introduce the
+        mutation.
+    :vartype start_time: float
+    :ivar single_site_id: The string ID of the single site at which to
+        introduce the mutation, see :meth:`Contig.add_single_site`.
+    :vartype single_site_id: str
+    :ivar population: The name of the population in which to introduce
+        the mutation.
+    :vartype population: str
     """
 
     time = attr.ib(type=float)
     single_site_id = attr.ib(type=str)
-    population_id = attr.ib(type=int)
+    population = attr.ib(type=str)
 
     def __attrs_post_init__(self):
         validate_time(self.time)
@@ -83,20 +86,38 @@ def __attrs_post_init__(self):
 @attr.s(kw_only=True)
 class ChangeMutationFitness(ExtendedEvent):
     """
-    TODO: docstring.
-
     Change the selection and/or dominance coefficients for the given mutation
     type in the given population. The mutation this applies to need not exist
     in the population for all generations between start_time and end_time.
     E.g. if the mutation could be transferred to the population via migration.
+    This requires that a `DrawMutation` event with the same `single_site_id` is
+    also supplied to the simulation engine.
 
     See SLiM documentation for registerFitnessCallback().
+
+    :ivar start_time: The number of generations in the past at which the change
+        in mutation fitness begins.
+    :vartype start_time: float
+    :ivar end_time: The number of generations in the past at which the change
+        in mutation fitness ends.
+    :vartype end_time: float
+    :ivar single_site_id: The string ID of the single site at which to change
+        the mutation fitness, see :meth:`Contig.add_single_site`.
+    :vartype single_site_id: str
+    :ivar population: The name of the population in which to change
+        the mutation fitness. If `None` (the default), then all populations are
+        affected.
+    :vartype population: str
+    :ivar selection_coeff: The new selection coefficient of the mutation.
+    :vartype selection_coeff: float
+    :ivar dominance_coeff: The new dominance coefficient of the mutation.
+    :vartype dominance_coeff: float
     """
 
     start_time = attr.ib(type=float)
     end_time = attr.ib(type=float)
     single_site_id = attr.ib(type=str)
-    population_id = attr.ib(type=int)
+    population = attr.ib(type=str, default=None)
     selection_coeff = attr.ib(type=float)
     dominance_coeff = attr.ib(type=float)
 
@@ -107,21 +128,35 @@ def __attrs_post_init__(self):
 @attr.s(kw_only=True)
 class ConditionOnAlleleFrequency(ExtendedEvent):
     """
-    TODO: docstring.
-
     Condition on the allele frequency of a drawn mutation with the given
     mutation type in the given population. The mutation need not have been
-    drawn in the population being conditioned upon.
-
-    This uses a save/restore mechanism in the SLiM code to do rejection
-    sampling on a simulation (when the condition is not met) without throwing
-    away the entire simulation.
+    drawn in the population being conditioned upon.  This requires that a
+    `DrawMutation` event for the same single site is also supplied to the
+    simulation engine.
+
+    :ivar start_time: The number of generations in the past at which to start
+        checking the allele frequency condition.
+    :vartype start_time: float
+    :ivar end_time: The number of generations in the past at which to stop
+        checking the allele frequency condition.
+    :vartype end_time: float
+    :ivar single_site_id: The string ID of the single site at which to check
+        the allele frequency condition, see :meth:`Contig.add_single_site`.
+    :vartype single_site_id: str
+    :ivar population: The name of the population in which to check
+        the allele frequency condition.
+    :vartype population: str
+    :ivar op: The comparison operator used for the allele frequency condition,
+        one of `(">", ">=", "<", "<=")`.
+    :vartype op: str
+    :ivar allele_frequency: The allele frequency to compare against.
+    :vartype op: float
     """
 
     start_time = attr.ib(type=float)
     end_time = attr.ib(type=float)
     single_site_id = attr.ib(type=str)
-    population_id = attr.ib(type=int)
+    population = attr.ib(type=str, default=None)
     op = attr.ib(type=str, default=None)
     allele_frequency = attr.ib(type=float, default=None)
 
@@ -151,3 +186,171 @@ def __attrs_post_init__(self):
     @classmethod
     def op_id(cls, op):
         return cls.op_types.index(op)
+
+
+def selective_sweep(
+    single_site_id,
+    population,
+    mutation_generation_ago,
+    selection_coeff,
+    dominance_coeff=0.5,
+    start_generation_ago=None,
+    end_generation_ago=None,
+    min_freq_at_start=None,
+    min_freq_at_end=None,
+    globally_adaptive=True,
+):
+    """
+    Creates a list of extended events corresponding to a selective sweep at a
+    single site (see :meth:`Contig.add_single_site`) that may subsequently be
+    passed to the simulation engine. The sequence of events is:
+
+        (1) In generation `mutation_generation_ago`, a **neutral** mutation is
+            introduced into `population`.
+
+        (2) In generation `start_generation_ago`, the mutation becomes
+            beneficial (whether globally or locally is controlled by the
+            `globally_adaptive` option). If the frequency of the mutation in
+            `population` is below `min_freq_at_start` in this generation, the
+            simulation will restart from `mutation_generation_ago`.
+
+        (3) In the generation after `end_generation_ago`, the mutation reverts to
+            neutrality in `population`. If the frequency of the mutation in
+            `population` is below `min_freq_at_end` in this generation, the
+            simulation will restart from `mutation_generation_ago`.
+
+    During the sweep, the fitnesses of individuals that are homozygous or
+    heterozygous for the beneficial mutation are `1 + selection_coeff` and  `1
+    + selection_coeff * dominance_coeff`, respectively. If `globally_adaptive`
+    is `True` (the default), then all individuals carrying the mutation are
+    impacted regardless of the population they belong to. Otherwise, the
+    mutation only influences fitness within `population`.
+
+    If the mutation is lost in `population` at any point between
+    `mutation_generation_ago` and `end_generation_ago`, then the simulation
+    will restart from `mutation_generation_ago`. The effect of "restarting"
+    (aka "rejection sampling") is that the resulting simulations are
+    conditioned on the mutation being extant after time
+    `mutation_generation_ago` with frequency above `min_freq_at_start` at
+    `start_generation_ago` and `min_freq_at_end` at `end_generation_ago`.
+
+    Because the simulation will restart until the allele frequency conditions
+    are satisfied, care must be taken to ensure that these are reasonable. For
+    example, if `start_generation_ago` is very close to
+    `mutation_generation_ago` and `min_freq_at_start` is substantially greater
+    than zero, the vast majority of proposed allele frequency trajectories will
+    be rejected and the simulation may take an infeasibly long time to
+    complete.
+
+    :param str single_site_id: The string ID of the single site at which to introduce
+        the mutation, see :meth:`Contig.add_single_site`.
+    :param int population: The name of the population in which the mutation is
+        introduced and the sweep occurs.
+    :param float mutation_generation_ago: The number of generations in the past at which
+        to introduce the mutation.
+    :param float selection_coeff: The selection coefficient of the beneficial
+        mutation during the sweep.  Must be greater than or equal to zero.
+    :param float dominance_coeff: The dominance coefficient of the beneficial
+        mutation during the sweep.  Defaults to `0.5`.
+    :param float start_generation_ago: The number of generations in the past at which the
+        sweep starts and the mutation switches from neutral to beneficial. Defaults to
+        `mutation_generation_ago` (i.e. a "hard" sweep).
+    :param float end_generation_ago: The number of generations in the past at which the
+        sweep ends and the mutation switches from beneficial to neutral. Defaults to
+        `0` (i.e. the end of the simulation).
+    :param float min_freq_at_start: The minimum allowed frequency for the
+        mutation in `population` at the start of the sweep (optional).
+    :param float min_freq_at_end: The minimum allowed frequency for the
+        mutation in `population` at the end of the sweep (optional).
+    :param bool globally_adaptive: If `True` (default), the mutation is beneficial in all
+        populations. If `False`, the mutation is beneficial only in the population of
+        origin.
+    """
+
+    if start_generation_ago is None:
+        start_generation_ago = mutation_generation_ago
+    if end_generation_ago is None:
+        end_generation_ago = 0
+
+    validate_time_range(mutation_generation_ago, start_generation_ago)
+    validate_time_range(start_generation_ago, end_generation_ago)
+    if selection_coeff < 0.0:
+        raise ValueError("Selection coefficient must be non-negative.")
+
+    extended_events = [
+        DrawMutation(
+            single_site_id=single_site_id,
+            time=mutation_generation_ago,
+            population=population,
+        ),
+        ConditionOnAlleleFrequency(
+            single_site_id=single_site_id,
+            start_time=GenerationAfter(mutation_generation_ago),
+            end_time=end_generation_ago,
+            population=population,
+            allele_frequency=0.0,
+            op=">",
+        ),
+    ]
+
+    # Sweep period is [start_generation_ago, end_generation_ago] (endpoints inclusive).
+    sweep_restricted_to = None if globally_adaptive else population
+    extended_events.extend(
+        [
+            ChangeMutationFitness(
+                single_site_id=single_site_id,
+                start_time=start_generation_ago,
+                end_time=end_generation_ago,
+                population=sweep_restricted_to,
+                selection_coeff=selection_coeff,
+                dominance_coeff=dominance_coeff,
+            ),
+        ]
+    )
+
+    # Only check sweep AF conditions in start/end generations, so that allele
+    # trajectories are otherwise unconstrained.
+    if min_freq_at_start is not None:
+        if not 0 < min_freq_at_start <= 1:
+            raise ValueError(
+                "If specified, the minimum allele frequency at the start of the sweep "
+                "must be in (0, 1]."
+            )
+        if start_generation_ago == mutation_generation_ago:
+            raise ValueError(
+                "If the start time of the sweep coincides with the introduction "
+                "of the mutation, then `min_freq_at_start` cannot be set."
+            )
+        extended_events.extend(
+            [
+                ConditionOnAlleleFrequency(
+                    single_site_id=single_site_id,
+                    start_time=start_generation_ago,
+                    end_time=start_generation_ago,
+                    population=population,
+                    allele_frequency=min_freq_at_start,
+                    op=">=",
+                ),
+            ]
+        )
+
+    if min_freq_at_end is not None:
+        if not 0 < min_freq_at_end <= 1:
+            raise ValueError(
+                "If specified, the minimum allele frequency at the end of the sweep "
+                "must be in (0, 1]."
+            )
+        extended_events.extend(
+            [
+                ConditionOnAlleleFrequency(
+                    single_site_id=single_site_id,
+                    start_time=end_generation_ago,
+                    end_time=end_generation_ago,
+                    population=population,
+                    allele_frequency=min_freq_at_end,
+                    op=">=",
+                ),
+            ]
+        )
+
+    return extended_events
diff --git a/stdpopsim/genomes.py b/stdpopsim/genomes.py
index fc14a1708..04543a8ee 100644
--- a/stdpopsim/genomes.py
+++ b/stdpopsim/genomes.py
@@ -484,52 +484,38 @@ def add_single_site(
         self,
         id,
         coordinate,
-        selection_coeff=0.0,
-        dominance_coeff=0.5,
         description=None,
         long_description=None,
     ):
         """
-        Adds a model of a single site mutation at the provided coordinate.
+        Adds a single site mutation class at the provided coordinate.
         The string label of the single site may be referenced by extended
         events passed to the simulation engine.
 
-        For instance, to simulate a selective sweep,
+        For instance, to simulate a selective sweep in population `pop_0`
+        starting at 1000 generations in the past:
 
         .. code-block:: python
 
             contig.add_single_site(
                 id="hard_sweep",
-                selection_coeff=0.1,
                 coordinate=mutation_coordinate,
             )
-            extended_events = [
-                stdpopsim.ext.DrawMutation(
-                    time=mutation_time,
-                    single_site_id="hard_sweep",
-                    population_id=0,
-                ),
-                stdpopsim.ext.ConditionOnAlleleFrequency(
-                    start_time=stdpopsim.ext.GenerationAfter(mutation_time),
-                    end_time=0,
-                    single_site_id="hard_sweep",
-                    population_id=0,
-                    op=">",
-                    allele_frequency=0.0,
-                ),
-            ]
+            extended_events = ext.selective_sweep(
+                mutation_generation_ago=1000,
+                selection_coeff=0.1,
+                population="pop_0",
+            )
             engine = stdpopsim.get_engine("slim")
             ts_sweep = engine.simulate(
                 ...,
                 extended_events=extended_events,
             )
 
+        See :func:`ext.selective_sweep` for more details on the sweep model.
+
         :param str id: A unique identifier for the single site.
         :param int coordinate: The coordinate of the site on the contig.
-        :param float selection_coeff: The selection coefficient of the mutation,
-            see :class:`MutationType` for details.
-        :param float dominance_coeff: The dominance coefficient of the mutation,
-            see :class:`MutationType` for details.
         :param str description: A short description of the single site mutation
             model as it would be used in written text, e.g., "Strong selective
             sweep".
@@ -542,10 +528,11 @@ def add_single_site(
 
         mut_type = stdpopsim.MutationType(
             distribution_type="f",
-            dominance_coeff=dominance_coeff,
-            distribution_args=[selection_coeff],
+            dominance_coeff=0.5,
+            distribution_args=[0.0],
             convert_to_substitution=False,
         )
+
         dfe = stdpopsim.DFE(
             id=id,
             mutation_types=[mut_type],
diff --git a/stdpopsim/slim_engine.py b/stdpopsim/slim_engine.py
index d529822ac..bd173eb84 100644
--- a/stdpopsim/slim_engine.py
+++ b/stdpopsim/slim_engine.py
@@ -51,6 +51,7 @@
 import textwrap
 import logging
 import warnings
+import collections
 
 import stdpopsim
 import numpy as np
@@ -210,11 +211,12 @@
      *
      * 1. There may be additional late events for the tick in which
      * restore() was called, and they are still scheduled to run.
-     * So we deactivate all script blocks to avoid unexpected problems.
-     * They will be automatically reactivated at the start of the next
-     * tick (see SLiM manual section 23.10).
+     * So we deactivate all script blocks in the "late" cycle to avoid
+     * unexpected problems. They will be automatically reactivated at the
+     * start of the next tick.
      */
-    community.allScriptBlocks.active = F;
+    sb = community.allScriptBlocks;
+    sb[sb.type == "late"].active = F;
 
     /*
      * 2. The late events below were run in the save() tick,
@@ -468,21 +470,36 @@
                     g_start+" > g_end="+g_end);
             }
 
-            community.registerLateEvent(NULL,
-                "{dbg('s="+selection_coeff+", h="+dominance_coeff+
-                " for m"+mut_type+" in p"+pop_id+"');}",
-                g_start, g_start);
-            community.registerLateEvent(NULL,
-                "{dbg('s, h defaults for m"+mut_type+" in p"+pop_id+"');}",
-                g_end, g_end);
             /* We explicitly format() here to prevent integral-valued floats
              * from getting converted to integers during string interpolation
              * (this triggers a type error when the fitness callback runs). */
             f_hom = format("%e", 1 + selection_coeff);
             f_het = format("%e", 1 + selection_coeff * dominance_coeff);
-            sim.registerMutationEffectCallback(NULL,
-                "{if (homozygous) return "+f_hom+"; else return "+f_het+";}",
-                mut_type, pop_id, g_start, g_end);
+
+            /* "All populations" is encoded by a negative value of pop_id. */
+            if (pop_id < 0) {
+                community.registerLateEvent(NULL,
+                    "{dbg('s="+selection_coeff+", h="+dominance_coeff+
+                    " for m"+mut_type+" globally');}",
+                    g_start, g_start);
+                community.registerLateEvent(NULL,
+                    "{dbg('s, h defaults for m"+mut_type+" globally');}",
+                    g_end, g_end);
+                sim.registerMutationEffectCallback(NULL,
+                    "{if (homozygous) return "+f_hom+"; else return "+f_het+";}",
+                    mut_type, NULL, g_start, g_end);
+            } else {
+                community.registerLateEvent(NULL,
+                    "{dbg('s="+selection_coeff+", h="+dominance_coeff+
+                    " for m"+mut_type+" in p"+pop_id+"');}",
+                    g_start, g_start);
+                community.registerLateEvent(NULL,
+                    "{dbg('s, h defaults for m"+mut_type+" in p"+pop_id+"');}",
+                    g_end, g_end);
+                sim.registerMutationEffectCallback(NULL,
+                    "{if (homozygous) return "+f_hom+"; else return "+f_het+";}",
+                    mut_type, pop_id, g_start, g_end);
+            }
         }
     }
 
@@ -530,6 +547,7 @@
 1 early () {
     defineConstant("log", community.createLogFile("$logfile", logInterval=NULL));
     log.addTick();
+    log.precision = 12;
     for (pop in sim.subpopulations.id) {
         log.addMeanSDColumns(
             "fitness_p" + pop,
@@ -539,7 +557,7 @@
 }
 
 1: early() {
-    if (community.tick % $loginterval == 1)
+    if ((community.tick - 1) % $loginterval == 0)
         log.logRow();
 }
 """
@@ -945,13 +963,17 @@ def fix_time(event):
     condition_on_allele_frequency = []
     op_id = stdpopsim.ext.ConditionOnAlleleFrequency.op_id
     slim_mutation_ids = _dfe_to_mtypes(contig)
-    drawn_single_site_ids = set()
+    drawn_single_site_ids = collections.defaultdict(int)
     referenced_single_site_ids = set()
+    population_name_to_index = {
+        x.name: i for i, x in enumerate(demographic_model.populations)
+    }
     for ee in extended_events:
+        cls_name = ee.__class__.__name__
         mutation_type_id = None
+        population_id = None
         coordinate = None
         if hasattr(ee, "single_site_id"):
-            cls_name = ee.__class__.__name__
             dfe_index = [
                 i for i, x in enumerate(contig.dfe_list) if x.id == ee.single_site_id
             ]
@@ -1001,38 +1023,54 @@ def fix_time(event):
             # Now that GenerationAfter times have been accounted for, we can
             # properly catch invalid start/end times.
             stdpopsim.ext.validate_time_range(ee.start_time, ee.end_time)
+        if hasattr(ee, "population"):
+            # Convert population name to integer index. "-1" is used to encode
+            # all populations (currently only valid for ChangeMutationFitness).
+            if ee.population is None:
+                population_id = -1
+            elif ee.population not in population_name_to_index.keys():
+                raise ValueError(
+                    f"The population {ee.population} referenced by "
+                    f"{cls_name} is not in demographic model "
+                    f"{demographic_model.id}, with defined populations "
+                    f"{population_name_to_index.keys()}."
+                )
+            else:
+                population_id = population_name_to_index[ee.population]
 
         # Append attributes to lists per event type
         if isinstance(ee, stdpopsim.ext.DrawMutation):
+            assert population_id >= 0
             drawn_mutations.append(
                 (
                     ee.time * demographic_model.generation_time,
                     mutation_type_id,
-                    ee.population_id,
+                    population_id,
                     coordinate,
                     0,  # flag to save state in mutation generation
                 )
             )
-            drawn_single_site_ids.add(ee.single_site_id)
+            drawn_single_site_ids[ee.single_site_id] += 1
         elif isinstance(ee, stdpopsim.ext.ChangeMutationFitness):
             fitness_callbacks.append(
                 (
                     ee.start_time * demographic_model.generation_time,
                     ee.end_time * demographic_model.generation_time,
                     mutation_type_id,
-                    ee.population_id,
+                    population_id,
                     ee.selection_coeff,
                     ee.dominance_coeff,
                 )
             )
             referenced_single_site_ids.add(ee.single_site_id)
         elif isinstance(ee, stdpopsim.ext.ConditionOnAlleleFrequency):
+            assert population_id >= 0
             condition_on_allele_frequency.append(
                 (
                     ee.start_time * demographic_model.generation_time,
                     ee.end_time * demographic_model.generation_time,
                     mutation_type_id,
-                    ee.population_id,
+                    population_id,
                     op_id(ee.op),
                     ee.allele_frequency,
                 )
@@ -1041,11 +1079,20 @@ def fix_time(event):
         else:
             raise ValueError(f"Unknown extended event type {type(ee)}")
 
+    # Check that there exists only one drawn mutation per single site.
+    for single_site_id, n_drawn_mutations in drawn_single_site_ids.items():
+        if n_drawn_mutations > 1:
+            raise ValueError(
+                f"The single site {single_site_id} has {n_drawn_mutations} "
+                f"mutations, but a maximum of one mutation is allowed per "
+                f"single site."
+            )
+
     # Check that drawn mutations exist for extended events that need them.
     for single_site_id in referenced_single_site_ids:
-        if single_site_id not in drawn_single_site_ids:
+        if single_site_id not in drawn_single_site_ids.keys():
             raise ValueError(
-                f"An extended event requires a mutation at single site id "
+                f"An extended event requires a mutation at single site "
                 f"{single_site_id}, but no mutation is drawn at this site."
             )
 
@@ -1420,6 +1467,9 @@ def simulate(
 
         :param seed: The seed for the random number generator.
         :type seed: int
+        :param extended_events: A list of :class:`ext.ExtendedEvents` to be
+            passed to SLiM, e.g. produced by :func:`ext.selective_sweep()`.
+        :type extended_events: list
         :param slim_path: The full path to the slim executable, or the name of
             a command in the current PATH.
         :type slim_path: str
@@ -1519,7 +1569,7 @@ def script_file_f():
         ts = _add_dfes_to_metadata(ts, contig)
         if _recap_and_rescale:
             ts = self._recap_and_rescale(
-                ts, seed, recap_epoch, contig, slim_scaling_factor
+                ts, seed, recap_epoch, contig, slim_scaling_factor, extended_events
             )
 
         if contig.inclusion_mask is not None:
@@ -1598,7 +1648,9 @@ def _simplify_remembered(self, ts):
         )
         return ts.simplify(samples=list(nodes), filter_populations=False)
 
-    def _recap_and_rescale(self, ts, seed, recap_epoch, contig, slim_scaling_factor):
+    def _recap_and_rescale(
+        self, ts, seed, recap_epoch, contig, slim_scaling_factor, extended_events=None
+    ):
         """
         Apply post-SLiM transformations to ``ts``. This rescales node times,
         does recapitation, simplification, and adds neutral mutations.
@@ -1610,6 +1662,17 @@ def _recap_and_rescale(self, ts, seed, recap_epoch, contig, slim_scaling_factor)
             table.time *= slim_scaling_factor
         metadata = tables.metadata
         metadata["SLiM"]["tick"] *= slim_scaling_factor
+        # If a single site DFE is referenced by an extended event, mark the DFE
+        # as non-neutral so that mutations aren't simulated over it.
+        if extended_events is not None:
+            single_sites = set()
+            for ee in extended_events:
+                assert hasattr(ee, "single_site_id")
+                single_sites.add(ee.single_site_id)
+            for dfe in metadata["stdpopsim"]["DFEs"]:
+                if dfe["id"] in single_sites:
+                    assert len(dfe["mutation_types"]) == 1
+                    dfe["mutation_types"][0]["is_neutral"] = False
         # Finding what slim id to use in recap DFE
         max_id = -1
         for dfe in metadata["stdpopsim"]["DFEs"]:
@@ -1762,7 +1825,9 @@ def recap_and_rescale(
                 slim_rate_map,
             )
 
-        ts = self._recap_and_rescale(ts, seed, recap_epoch, contig, slim_scaling_factor)
+        ts = self._recap_and_rescale(
+            ts, seed, recap_epoch, contig, slim_scaling_factor, extended_events
+        )
         return ts
 
 


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 98aeb688d..2fdd1d6a6 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -17,7 +17,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-18.04, macos-10.15, windows-latest]
+        os: [ubuntu-20.04, macos-10.15, windows-latest]
         python: [3.7, "3.10"]
     env:
       CONDA_ENV_NAME: stdpopsim
diff --git a/tests/test_slim_engine.py b/tests/test_slim_engine.py
index c879bb19b..5814c7e3e 100644
--- a/tests/test_slim_engine.py
+++ b/tests/test_slim_engine.py
@@ -9,6 +9,7 @@
 from unittest import mock
 import numpy as np
 from numpy.testing import assert_array_equal
+import collections
 
 import pytest
 import tskit
@@ -1363,7 +1364,7 @@ def test_draw_mutation(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id=self.mut_id,
-                population_id=0,
+                population="pop_0",
             ),
         ]
         engine = stdpopsim.get_engine("slim")
@@ -1375,20 +1376,20 @@ def test_draw_mutation(self):
             dry_run=True,
         )
 
-    def test_draw_multiple_mutations(self):
+    def test_draw_mutations_at_different_sites(self):
         contig = get_test_contig()
         contig.add_single_site(id="recent", coordinate=100)
         contig.add_single_site(id="older", coordinate=101)
         extended_events = [
             stdpopsim.ext.DrawMutation(
-                time=int(self.T_mut / 2),
+                time=self.T_mut // 2,
                 single_site_id="recent",
-                population_id=0,
+                population="pop_0",
             ),
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id="older",
-                population_id=0,
+                population="pop_0",
             ),
         ]
         engine = stdpopsim.get_engine("slim")
@@ -1400,6 +1401,31 @@ def test_draw_multiple_mutations(self):
             dry_run=True,
         )
 
+    def test_draw_multiple_mutations_at_same_site(self):
+        contig = get_test_contig()
+        contig.add_single_site(id="mutant", coordinate=100)
+        extended_events = [
+            stdpopsim.ext.DrawMutation(
+                time=self.T_mut // 2,
+                single_site_id="mutant",
+                population="pop_0",
+            ),
+            stdpopsim.ext.DrawMutation(
+                time=self.T_mut,
+                single_site_id="mutant",
+                population="pop_0",
+            ),
+        ]
+        engine = stdpopsim.get_engine("slim")
+        with pytest.raises(ValueError, match="maximum of one mutation is allowed"):
+            engine.simulate(
+                demographic_model=self.model,
+                contig=contig,
+                samples=self.samples,
+                extended_events=extended_events,
+                dry_run=True,
+            )
+
     def test_invalid_single_site_id(self):
         engine = stdpopsim.get_engine("slim")
         for single_site_id in ["deleterious", "sweep"]:
@@ -1407,7 +1433,7 @@ def test_invalid_single_site_id(self):
                 stdpopsim.ext.DrawMutation(
                     time=self.T_mut,
                     single_site_id=single_site_id,
-                    population_id=0,
+                    population="pop_0",
                 ),
             ]
             with pytest.raises(ValueError, match="must exist and be uniquely labelled"):
@@ -1424,7 +1450,7 @@ def test_no_mutation_types_defined(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id=self.mut_id,
-                population_id=0,
+                population="pop_0",
             ),
         ]
         contig = get_test_contig()
@@ -1445,7 +1471,7 @@ def test_multiple_mutation_types_defined(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
             ),
         ]
         contig = get_test_contig()
@@ -1481,7 +1507,7 @@ def test_fitness_distribution_not_fixed(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
             ),
         ]
         contig = get_test_contig()
@@ -1508,7 +1534,7 @@ def test_multiple_intervals(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
             ),
         ]
         contig = get_test_contig()
@@ -1544,7 +1570,7 @@ def test_interval_too_large(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
             ),
         ]
         contig = get_test_contig()
@@ -1580,7 +1606,7 @@ def test_duplicate_mutation_ids(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
             ),
         ]
         contig = get_test_contig()
@@ -1601,7 +1627,7 @@ def test_mutation_has_no_interval(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
             ),
         ]
         contig = get_test_contig()
@@ -1623,14 +1649,14 @@ def test_bad_time(self):
                 stdpopsim.ext.DrawMutation(
                     time=time,
                     single_site_id="irrelevant",
-                    population_id=0,
+                    population="pop_0",
                 )
         for time in (0, -1):
             with pytest.raises(ValueError):
                 stdpopsim.ext.DrawMutation(
                     time=stdpopsim.ext.GenerationAfter(time),
                     single_site_id="irrelevant",
-                    population_id=0,
+                    population="pop_0",
                 )
 
 
@@ -1646,13 +1672,13 @@ def test_drawn_mutation_not_lost(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
             ),
             stdpopsim.ext.ConditionOnAlleleFrequency(
                 start_time=0,
                 end_time=0,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
                 op=">",
                 allele_frequency=0,
             ),
@@ -1678,13 +1704,13 @@ def test_drawn_mutation_is_lost(self):
             stdpopsim.ext.DrawMutation(
                 time=self.T_mut,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
             ),
             stdpopsim.ext.ConditionOnAlleleFrequency(
                 start_time=0,
                 end_time=0,
                 single_site_id="test",
-                population_id=0,
+                population="pop_0",
                 op="<=",
                 allele_frequency=0,
             ),
@@ -1711,14 +1737,14 @@ def test_drawn_mutation_meets_AF_threshold(self):
                 stdpopsim.ext.DrawMutation(
                     time=self.T_mut,
                     single_site_id="test",
-                    population_id=0,
+                    population="pop_0",
                 ),
                 # Condition on desired AF at end of simulation.
                 stdpopsim.ext.ConditionOnAlleleFrequency(
                     start_time=0,
                     end_time=0,
                     single_site_id="test",
-                    population_id=0,
+                    population="pop_0",
                     op=">=",
                     allele_frequency=af_threshold,
                 ),
@@ -1759,7 +1785,7 @@ def test_bad_AF_conditioning_parameters(self):
                     start_time=0,
                     end_time=0,
                     single_site_id="irrelevant",
-                    population_id=0,
+                    population="pop_0",
                     op=op,
                     allele_frequency=af,
                 )
@@ -1771,7 +1797,7 @@ def test_bad_times(self):
                     start_time=start_time,
                     end_time=end_time,
                     single_site_id="irrelevant",
-                    population_id=0,
+                    population="pop_0",
                     op=">",
                     allele_frequency=0,
                 )
@@ -1794,13 +1820,13 @@ def test_bad_GenerationAfter_times(self):
                     stdpopsim.ext.DrawMutation(
                         time=self.T_mut,
                         single_site_id=self.mut_id,
-                        population_id=0,
+                        population="pop_0",
                     ),
                     stdpopsim.ext.ConditionOnAlleleFrequency(
                         start_time=stdpopsim.ext.GenerationAfter(start_time),
                         end_time=end_time,
                         single_site_id=self.mut_id,
-                        population_id=0,
+                        population="pop_0",
                         op=">",
                         allele_frequency=0,
                     ),
@@ -1828,7 +1854,7 @@ def test_no_drawn_mutation(self):
                 start_time=stdpopsim.ext.GenerationAfter(self.T_mut),
                 end_time=0,
                 single_site_id=self.mut_id,
-                population_id=0,
+                population="pop_0",
                 op=">",
                 allele_frequency=0,
             ),
@@ -1844,51 +1870,6 @@ def test_no_drawn_mutation(self):
             )
 
 
-@pytest.mark.skipif(IS_WINDOWS, reason="SLiM not available on windows")
-class TestFixedSelectionCoefficient(PiecewiseConstantSizeMixin):
-    @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
-    def test_drawn_mutation_has_correct_selection_coeff(self):
-        engine = stdpopsim.get_engine("slim")
-        contig = get_test_contig()
-        mut_id = "sweep"
-        contig.add_single_site(
-            id=mut_id,
-            coordinate=100,
-            selection_coeff=0.1,
-            dominance_coeff=0.5,
-        )
-        extended_events = [
-            stdpopsim.ext.DrawMutation(
-                time=self.T_mut,
-                single_site_id=mut_id,
-                population_id=0,
-            ),
-            stdpopsim.ext.ConditionOnAlleleFrequency(
-                start_time=0,
-                end_time=0,
-                single_site_id=mut_id,
-                population_id=0,
-                op=">",
-                allele_frequency=0,
-            ),
-        ]
-        scaling_factor = 10
-        ts = engine.simulate(
-            demographic_model=self.model,
-            contig=contig,
-            samples=self.samples,
-            extended_events=extended_events,
-            slim_scaling_factor=scaling_factor,
-            slim_burn_in=0.1,
-            _recap_and_rescale=False,
-        )
-        assert ts.num_mutations == 1
-        mut = next(ts.mutations())
-        mutation_list = mut.metadata["mutation_list"]
-        assert len(mutation_list) == 1
-        assert mutation_list[0]["selection_coeff"] == scaling_factor * 0.1
-
-
 @pytest.mark.skipif(IS_WINDOWS, reason="SLiM not available on windows")
 class TestChangeMutationFitness(PiecewiseConstantSizeMixin):
     # Testing stdpopsim.ext.ChangeMutationFitness is challenging, because
@@ -1906,13 +1887,13 @@ def test_positive_mutation_meets_AF_threshold(self):
                 stdpopsim.ext.DrawMutation(
                     time=self.T_mut,
                     single_site_id=self.mut_id,
-                    population_id=0,
+                    population="pop_0",
                 ),
                 stdpopsim.ext.ChangeMutationFitness(
                     start_time=stdpopsim.ext.GenerationAfter(self.T_mut),
                     end_time=0,
                     single_site_id=self.mut_id,
-                    population_id=0,
+                    population="pop_0",
                     selection_coeff=0.1,
                     dominance_coeff=0.5,
                 ),
@@ -1922,7 +1903,7 @@ def test_positive_mutation_meets_AF_threshold(self):
                     start_time=stdpopsim.ext.GenerationAfter(self.T_mut),
                     end_time=0,
                     single_site_id=self.mut_id,
-                    population_id=0,
+                    population="pop_0",
                     op=">",
                     allele_frequency=0,
                 ),
@@ -1931,7 +1912,7 @@ def test_positive_mutation_meets_AF_threshold(self):
                     start_time=0,
                     end_time=0,
                     single_site_id=self.mut_id,
-                    population_id=0,
+                    population="pop_0",
                     op=">=",
                     allele_frequency=af_threshold,
                 ),
@@ -1949,13 +1930,38 @@ def test_positive_mutation_meets_AF_threshold(self):
             assert ts.num_mutations == 1
             assert self.allele_frequency(ts) >= af_threshold
 
+    @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
+    def test_referenced_single_site_is_nonneutral(self):
+        engine = stdpopsim.get_engine("slim")
+        contig = get_test_contig()
+        contig.add_single_site("one", coordinate=100)
+        contig.add_single_site("two", coordinate=101)
+        extended_events = [
+            stdpopsim.ext.DrawMutation(
+                time=self.T_mut,
+                single_site_id="one",
+                population="pop_0",
+            ),
+        ]
+        ts = engine.simulate(
+            demographic_model=self.model,
+            contig=contig,
+            samples=self.samples,
+            extended_events=extended_events,
+        )
+        referenced_dfe = ts.metadata["stdpopsim"]["DFEs"][1]
+        assert referenced_dfe["id"] == "one"
+        assert referenced_dfe["mutation_types"][0]["is_neutral"] is False
+        unreferenced_dfe = ts.metadata["stdpopsim"]["DFEs"][2]
+        assert unreferenced_dfe["id"] == "two"
+        assert unreferenced_dfe["mutation_types"][0]["is_neutral"] is True
+
     def test_no_drawn_mutation(self):
         extended_events = [
             stdpopsim.ext.ChangeMutationFitness(
                 start_time=stdpopsim.ext.GenerationAfter(self.T_mut),
                 end_time=0,
                 single_site_id=self.mut_id,
-                population_id=0,
                 selection_coeff=0.1,
                 dominance_coeff=0.5,
             ),
@@ -1977,11 +1983,33 @@ def test_bad_times(self):
                     start_time=start_time,
                     end_time=end_time,
                     single_site_id="irrelevant",
-                    population_id=0,
+                    population="pop_0",
                     selection_coeff=0.1,
                     dominance_coeff=0.5,
                 )
 
+    def test_population_name(self):
+        engine = stdpopsim.get_engine("slim")
+        for pop in ["not present", 0]:
+            extended_events = [
+                stdpopsim.ext.ChangeMutationFitness(
+                    start_time=stdpopsim.ext.GenerationAfter(self.T_mut),
+                    end_time=0,
+                    single_site_id=self.mut_id,
+                    population=pop,
+                    selection_coeff=0.1,
+                    dominance_coeff=0.5,
+                ),
+            ]
+            with pytest.raises(ValueError, match="is not in demographic model"):
+                engine.simulate(
+                    demographic_model=self.model,
+                    contig=self.contig,
+                    samples=self.samples,
+                    extended_events=extended_events,
+                    dry_run=True,
+                )
+
     def test_bad_GenerationAfter_times(self):
         engine = stdpopsim.get_engine("slim")
         for start_time, end_time in [
@@ -2000,13 +2028,13 @@ def test_bad_GenerationAfter_times(self):
                     stdpopsim.ext.DrawMutation(
                         time=self.T_mut,
                         single_site_id=self.mut_id,
-                        population_id=0,
+                        population="pop_0",
                     ),
                     stdpopsim.ext.ChangeMutationFitness(
                         start_time=stdpopsim.ext.GenerationAfter(start_time),
                         end_time=end_time,
                         single_site_id=self.mut_id,
-                        population_id=0,
+                        population="pop_0",
                         selection_coeff=0.1,
                         dominance_coeff=0.5,
                     ),
@@ -2038,3 +2066,449 @@ def test_bad_extended_events(self):
                     extended_events=[bad_ee],
                     dry_run=True,
                 )
+
+
+@pytest.mark.skipif(IS_WINDOWS, reason="SLiM not available on windows")
+class TestSelectiveSweep(PiecewiseConstantSizeMixin):
+    @staticmethod
+    def _get_island_model(Ne=1000, migration_rate=0.01):
+        model = msprime.Demography()
+        model.add_population(initial_size=Ne, name="pop0")
+        model.add_population(initial_size=Ne, name="pop1")
+        model.set_migration_rate(source="pop0", dest="pop1", rate=migration_rate)
+        model.set_migration_rate(source="pop1", dest="pop0", rate=migration_rate)
+        return stdpopsim.DemographicModel(
+            id="🏝️",
+            description="island model",
+            long_description="for sweep tests",
+            model=model,
+            generation_time=1,
+        )
+
+    @staticmethod
+    def _fitness_per_generation(
+        logfile, start_generation_ago, end_generation_ago, pop=0
+    ):
+        # NB: Be careful of rounding error in log with weak selection
+        assert start_generation_ago >= end_generation_ago
+
+        # Pull out average fitness for every generation in a specified time
+        # interval from the fitness logfile. Columns are:
+        #   "tick", "fitness_p0_mean", "fitness_p0_sd"
+        # This test will break if the logfile format changes.
+        log = np.genfromtxt(logfile, delimiter=",", names=True)
+        tick = log["tick"].astype(int)
+
+        # Because of save-restore during rejection sampling, the log may
+        # contain multiple trajectories; keep the one that passed (last
+        # fitness value for a given generation).
+        accepted_trajectory = {
+            i: x for i, x in zip(tick, log["fitness_p" + str(pop) + "_mean"])
+        }
+        generation = np.array(sorted(accepted_trajectory))
+        mean_fitness = np.array([accepted_trajectory[i] for i in generation])
+        generation = np.max(generation) - generation
+
+        # Generations in which trajectory was rejected
+        rejected = [tick[i + 1] <= tick[i] for i in range(len(tick) - 1)] + [False]
+        rejections_by_generation = collections.Counter(np.max(tick) - tick[rejected])
+
+        # Check that logging interval hits sweep boundaries exactly.  Fitness
+        # is calculated in the early() block; so in the logfile, mean fitness
+        # refers to the fitness of the generation **previous** to the tick.
+        start_generation_ago = max(start_generation_ago - 1, 0)
+        end_generation_ago = max(end_generation_ago - 1, 0)
+        in_sweep = np.logical_and(
+            generation <= start_generation_ago,
+            generation >= end_generation_ago,
+        )
+        assert np.max(generation[in_sweep]) == start_generation_ago
+        assert np.min(generation[in_sweep]) == end_generation_ago
+        assert np.sum(in_sweep) == start_generation_ago - end_generation_ago + 1
+
+        return (
+            mean_fitness[in_sweep],
+            mean_fitness[~in_sweep],
+            rejections_by_generation,
+        )
+
+    @pytest.mark.usefixtures("tmp_path")
+    def test_sweep(self, tmp_path):
+        mutation_generation_ago = 1000
+        for start_generation_ago in [
+            mutation_generation_ago // 2,
+            mutation_generation_ago,
+        ]:
+            for end_generation_ago in [0, mutation_generation_ago // 4]:
+                logfile = tmp_path / "sweep.logfile"
+                engine = stdpopsim.get_engine("slim")
+                contig = get_test_contig()
+                locus_id = "sweep"
+                contig.add_single_site(
+                    id=locus_id,
+                    coordinate=100,
+                )
+                extended_events = stdpopsim.ext.selective_sweep(
+                    single_site_id=locus_id,
+                    population="pop_0",
+                    mutation_generation_ago=mutation_generation_ago,
+                    start_generation_ago=start_generation_ago,
+                    end_generation_ago=end_generation_ago,
+                    selection_coeff=0.1,
+                )
+                engine.simulate(
+                    demographic_model=self.model,
+                    contig=contig,
+                    samples=self.samples,
+                    extended_events=extended_events,
+                    slim_burn_in=1,
+                    logfile=logfile,
+                    logfile_interval=1,
+                )
+                in_sweep, outside_sweep, _ = self._fitness_per_generation(
+                    logfile=logfile,
+                    start_generation_ago=start_generation_ago,
+                    end_generation_ago=end_generation_ago,
+                    pop=0,
+                )
+                assert np.all(in_sweep > 1)
+                assert np.all(outside_sweep == 1)
+
+    @pytest.mark.usefixtures("tmp_path")
+    def test_sweep_meets_min_freq_at_start(self, tmp_path):
+        mutation_generation_ago = 10
+        start_generation_ago = 8
+        s = 0.01
+        # condition on a difficult-to-reach frequency, i.e.
+        min_freq = 3 * (mutation_generation_ago - start_generation_ago) / (2 * self.N0)
+        logfile = tmp_path / "sweep_start_af.logfile"
+        engine = stdpopsim.get_engine("slim")
+        contig = get_test_contig()
+        locus_id = "sweep"
+        contig.add_single_site(
+            id=locus_id,
+            coordinate=100,
+        )
+        extended_events = stdpopsim.ext.selective_sweep(
+            single_site_id=locus_id,
+            population="pop_0",
+            mutation_generation_ago=mutation_generation_ago,
+            start_generation_ago=start_generation_ago,
+            min_freq_at_start=min_freq,
+            selection_coeff=s,
+            dominance_coeff=0.5,
+        )
+        while True:
+            engine.simulate(
+                demographic_model=self.model,
+                contig=contig,
+                samples=self.samples,
+                extended_events=extended_events,
+                slim_burn_in=1,
+                logfile=logfile,
+                logfile_interval=1,
+            )
+            in_sweep, outside_sweep, rejections = self._fitness_per_generation(
+                logfile=logfile,
+                start_generation_ago=start_generation_ago,
+                end_generation_ago=0,
+                pop=0,
+            )
+            # ensure that rejections are occuring in the generation of the AF
+            # condition
+            if start_generation_ago in rejections.keys():
+                break
+        assert np.all(outside_sweep == 1.0)
+        assert in_sweep[0] >= s * min_freq + 1
+
+    @pytest.mark.usefixtures("tmp_path")
+    def test_sweep_meets_min_freq_at_end(self, tmp_path):
+        mutation_generation_ago = 10
+        end_generation_ago = 8
+        s = 0.01
+        # condition on a difficult-to-reach frequency, i.e.
+        min_freq = 3 * (mutation_generation_ago - end_generation_ago) / (2 * self.N0)
+        logfile = tmp_path / "sweep_end_af.logfile"
+        engine = stdpopsim.get_engine("slim")
+        contig = get_test_contig()
+        locus_id = "sweep"
+        contig.add_single_site(
+            id=locus_id,
+            coordinate=100,
+        )
+        extended_events = stdpopsim.ext.selective_sweep(
+            single_site_id=locus_id,
+            population="pop_0",
+            mutation_generation_ago=mutation_generation_ago,
+            end_generation_ago=end_generation_ago,
+            min_freq_at_end=min_freq,
+            selection_coeff=s,
+            dominance_coeff=0.5,
+        )
+        while True:
+            engine.simulate(
+                demographic_model=self.model,
+                contig=contig,
+                samples=self.samples,
+                extended_events=extended_events,
+                slim_burn_in=1,
+                logfile=logfile,
+                logfile_interval=1,
+            )
+            in_sweep, outside_sweep, rejections = self._fitness_per_generation(
+                logfile=logfile,
+                start_generation_ago=mutation_generation_ago,
+                end_generation_ago=end_generation_ago,
+                pop=0,
+            )
+            # ensure that rejections are occuring in the generation of the AF
+            # condition
+            if end_generation_ago in rejections.keys():
+                break
+        assert np.all(outside_sweep == 1.0)
+        assert in_sweep[-1] >= s * min_freq + 1
+
+    def test_sweep_with_negative_selection_coeff(self):
+        with pytest.raises(ValueError, match="coefficient must be"):
+            stdpopsim.ext.selective_sweep(
+                single_site_id="irrelevant",
+                population="irrelevant",
+                mutation_generation_ago=1000,
+                selection_coeff=-0.1,
+            )
+
+    def test_sweep_with_bad_AF_conditions(self):
+        for start_freq, end_freq in zip([-0.1, 0.1], [0.1, -0.1]):
+            with pytest.raises(ValueError, match="of the sweep must be in"):
+                stdpopsim.ext.selective_sweep(
+                    single_site_id="irrelevant",
+                    population="irrelevant",
+                    mutation_generation_ago=1000,
+                    selection_coeff=0.1,
+                    start_generation_ago=100,
+                    min_freq_at_start=start_freq,
+                    min_freq_at_end=end_freq,
+                )
+        with pytest.raises(
+            ValueError, match="coincides with the introduction of the mutation"
+        ):
+            stdpopsim.ext.selective_sweep(
+                single_site_id="irrelevant",
+                population="irrelevant",
+                mutation_generation_ago=1000,
+                start_generation_ago=1000,
+                selection_coeff=0.1,
+                min_freq_at_start=0.1,
+            )
+
+    @pytest.mark.usefixtures("tmp_path")
+    def test_global_sweep(self, tmp_path):
+        mutation_generation_ago = 100
+        start_generation_ago = 50
+        end_generation_ago = 30
+        logfile = tmp_path / "sweep_not_restricted.logfile"
+        engine = stdpopsim.get_engine("slim")
+        model = self._get_island_model()
+        contig = get_test_contig()
+        locus_id = "sweep"
+        contig.add_single_site(
+            id=locus_id,
+            coordinate=100,
+        )
+        extended_events = stdpopsim.ext.selective_sweep(
+            single_site_id=locus_id,
+            population="pop1",
+            mutation_generation_ago=mutation_generation_ago,
+            start_generation_ago=start_generation_ago,
+            end_generation_ago=end_generation_ago,
+            selection_coeff=1000.0,  # so that mutation stays in non-source pop
+            globally_adaptive=True,
+        )
+        while True:
+            # The while loop is to avoid rare stochastic failure due to
+            # mutation never migrating into the non-source population
+            engine.simulate(
+                demographic_model=model,
+                contig=contig,
+                samples=self.samples,
+                extended_events=extended_events,
+                slim_burn_in=1,
+                logfile=logfile,
+                logfile_interval=1,
+            )
+            p0_in_sweep, p0_outside_sweep, _ = self._fitness_per_generation(
+                logfile=logfile,
+                start_generation_ago=start_generation_ago,
+                end_generation_ago=end_generation_ago,
+                pop=0,
+            )
+            # Mutation is in non-source pop at start of sweep
+            if p0_in_sweep[0] > 1:
+                break
+        p1_in_sweep, p1_outside_sweep, _ = self._fitness_per_generation(
+            logfile=logfile,
+            start_generation_ago=start_generation_ago,
+            end_generation_ago=end_generation_ago,
+            pop=1,
+        )
+        assert np.all(p0_in_sweep > 1)
+        assert np.all(p0_outside_sweep == 1)
+        assert np.all(p1_in_sweep > 1)
+        assert np.all(p1_outside_sweep == 1)
+
+    @pytest.mark.usefixtures("tmp_path")
+    def test_local_sweep(self, tmp_path):
+        mutation_generation_ago = 100
+        start_generation_ago = 50
+        end_generation_ago = 30
+        logfile = tmp_path / "sweep_restricted.logfile"
+        engine = stdpopsim.get_engine("slim")
+        model = self._get_island_model()
+        contig = get_test_contig()
+        locus_id = "sweep"
+        contig.add_single_site(
+            id=locus_id,
+            coordinate=100,
+        )
+        extended_events = stdpopsim.ext.selective_sweep(
+            single_site_id=locus_id,
+            population="pop1",
+            mutation_generation_ago=mutation_generation_ago,
+            start_generation_ago=start_generation_ago,
+            end_generation_ago=end_generation_ago,
+            selection_coeff=0.1,
+            globally_adaptive=False,
+        )
+        engine.simulate(
+            demographic_model=model,
+            contig=contig,
+            samples=self.samples,
+            extended_events=extended_events,
+            slim_burn_in=1,
+            logfile=logfile,
+            logfile_interval=1,
+        )
+        p0_in_sweep, p0_outside_sweep, _ = self._fitness_per_generation(
+            logfile=logfile,
+            start_generation_ago=start_generation_ago,
+            end_generation_ago=end_generation_ago,
+            pop=0,
+        )
+        p1_in_sweep, p1_outside_sweep, _ = self._fitness_per_generation(
+            logfile=logfile,
+            start_generation_ago=start_generation_ago,
+            end_generation_ago=end_generation_ago,
+            pop=1,
+        )
+        assert np.all(p0_in_sweep == 1)
+        assert np.all(p0_outside_sweep == 1)
+        assert np.all(p1_in_sweep > 1)
+        assert np.all(p1_outside_sweep == 1)
+
+    @pytest.mark.usefixtures("tmp_path")
+    def test_sweeps_at_multiple_sites(self, tmp_path):
+        logfile = tmp_path / "sweep_multi.logfile"
+        engine = stdpopsim.get_engine("slim")
+        contig = get_test_contig()
+        model = self._get_island_model()
+        mutation_generation_ago = 100
+        ids = ["sweep1", "sweep2"]
+        coords = [100, contig.length - 100]
+        start_generation_agos = [mutation_generation_ago, mutation_generation_ago // 2]
+        end_generation_agos = [0, 0]
+        pop_ids = ["pop0", "pop1"]
+        extended_events = []
+        for locus_id, coord, start_generation_ago, end_generation_ago, pop in zip(
+            ids, coords, start_generation_agos, end_generation_agos, pop_ids
+        ):
+            contig.add_single_site(
+                id=locus_id,
+                coordinate=coord,
+            )
+            extended_events += stdpopsim.ext.selective_sweep(
+                single_site_id=locus_id,
+                population=pop,
+                mutation_generation_ago=mutation_generation_ago,
+                start_generation_ago=start_generation_ago,
+                end_generation_ago=end_generation_ago,
+                selection_coeff=1.0,
+                globally_adaptive=False,
+            )
+        engine.simulate(
+            demographic_model=model,
+            contig=contig,
+            samples=self.samples,
+            extended_events=extended_events,
+            slim_burn_in=1,
+            logfile=logfile,
+            logfile_interval=1,
+        )
+        for i, _ in enumerate(pop_ids):
+            in_sweep, outside_sweep, _ = self._fitness_per_generation(
+                logfile=logfile,
+                start_generation_ago=start_generation_agos[i],
+                end_generation_ago=end_generation_agos[i],
+                pop=i,
+            )
+            assert np.all(in_sweep > 1)
+            assert np.all(outside_sweep == 1)
+
+    @pytest.mark.usefixtures("tmp_path")
+    def test_sweep_with_background_selection(self, tmp_path):
+        logfile = tmp_path / "sweep_bgs.logfile"
+        engine = stdpopsim.get_engine("slim")
+        model = self._get_island_model()
+        contig = get_test_contig()
+        contig.add_dfe(
+            intervals=np.array([[0, contig.length // 2]], dtype="int"),
+            DFE=stdpopsim.DFE(
+                id="BGS",
+                description="deleterious",
+                long_description="mutations",
+                mutation_types=[
+                    stdpopsim.MutationType(
+                        distribution_type="e",
+                        distribution_args=[-0.01],
+                    )
+                ],
+            ),
+        )
+        locus_id = "sweep"
+        mutation_generation_ago = 1000
+        contig.add_single_site(
+            id=locus_id,
+            coordinate=100,
+        )
+        extended_events = stdpopsim.ext.selective_sweep(
+            single_site_id=locus_id,
+            population="pop1",
+            mutation_generation_ago=mutation_generation_ago,
+            selection_coeff=100.0,
+            globally_adaptive=False,
+        )
+        engine.simulate(
+            demographic_model=model,
+            contig=contig,
+            samples=self.samples,
+            extended_events=extended_events,
+            slim_burn_in=1,
+            logfile=logfile,
+            logfile_interval=1,
+        )
+        p0_in_sweep, p0_outside_sweep, _ = self._fitness_per_generation(
+            logfile=logfile,
+            start_generation_ago=mutation_generation_ago,
+            end_generation_ago=0,
+            pop=0,
+        )
+        p1_in_sweep, p1_outside_sweep, _ = self._fitness_per_generation(
+            logfile=logfile,
+            start_generation_ago=mutation_generation_ago,
+            end_generation_ago=0,
+            pop=1,
+        )
+        # Population 0 = BGS only; Population 1 = BGS + sweep
+        assert np.all(p0_in_sweep <= 1)
+        assert np.all(p0_outside_sweep <= 1)
+        assert np.all(p1_in_sweep > 1)
+        assert np.all(p1_outside_sweep <= 1)