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

Issue:
Cannot `recap_and_rescale` tree sequence produced by saved SLiM script, because DFE metadata is missing
The "Outputting a SLiM script" section of the tutorial (a few subsections down from [here](https://popsim-consortium.github.io/stdpopsim-docs/latest/tutorial.html#an-example-simulation)) gives an example of saving a SLiM script, then loading it back into python and recapitating via `Engine.recap_and_rescale`.

But, this doesn't work:
```
ts = tskit.load("foo.trees")
ts = engine.recap_and_rescale(ts, demography, contig, samples, slim_scaling_factor=10)                                                                                          
# Traceback (most recent call last):
#   File "<stdin>", line 1, in <module>
#   File "/home/natep/projects/stdpopsim/stdpopsim/stdpopsim/slim_engine.py", line 1900, in recap_and_rescale
#     ts = self._recap_and_rescale(
#   File "/home/natep/projects/stdpopsim/stdpopsim/stdpopsim/slim_engine.py", line 1719, in _recap_and_rescale
#     for dfe in metadata["stdpopsim"]["DFEs"]:
# KeyError: 'stdpopsim'
```

The reason is because the DFE metadata is not created within the SLiM script -- it's retroactively added by `simulate` here:
https://github.com/popsim-consortium/stdpopsim/blob/e0900f1a41224b605faba6c9409966ac56b5c93f/stdpopsim/slim_engine.py#L1596

I think what we need to do here is add the same line to `Engine.recap_and_rescale` (not `Engine._recap_and_rescale` which is also used by `Engine.simulate`).


Discussion / hints:
Should probably also add a test that:
- Saves a SLiM script
- Uses sed to alter output path for trees file to a non-temporary location
- Runs SLiM
- Loads in trees and recapitates

Fix patch:
diff --git a/stdpopsim/slim_engine.py b/stdpopsim/slim_engine.py
index 58189b6c1..6de3adf99 100644
--- a/stdpopsim/slim_engine.py
+++ b/stdpopsim/slim_engine.py
@@ -1897,6 +1897,7 @@ def recap_and_rescale(
                 slim_rate_map,
             )
 
+        ts = _add_dfes_to_metadata(ts, contig)
         ts = self._recap_and_rescale(
             ts,
             seed,


Test patch:
diff --git a/tests/test_slim_engine.py b/tests/test_slim_engine.py
index 2db8025c7..537afac3e 100644
--- a/tests/test_slim_engine.py
+++ b/tests/test_slim_engine.py
@@ -10,6 +10,7 @@
 import numpy as np
 from numpy.testing import assert_array_equal
 import collections
+import re
 
 import pytest
 import tskit
@@ -266,6 +267,60 @@ def test_recap_and_rescale(self):
                 assert tables1.edges == tables2.edges
                 assert tables1.mutations == tables2.mutations
 
+    @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
+    @pytest.mark.usefixtures("tmp_path")
+    def test_recap_and_rescale_on_external_slim_run(self, tmp_path):
+        # _SLiMEngine.simulate() adds metadata after SLiM is run.  If we
+        # produce a script file, then run with SLiM externally, this metadata
+        # would be omitted. We want to check that engine.recap_and_rescale
+        # still runs in this case.
+        treefile = tmp_path / "foo.trees"
+        scriptfile = tmp_path / "slim.script"
+        engine = stdpopsim.get_engine("slim")
+        species = stdpopsim.get_species("HomSap")
+        contig = species.get_contig("chr1", length_multiplier=0.01)
+        model = stdpopsim.PiecewiseConstantSize(species.population_size)
+        samples = {"pop_0": 5}
+        seed = 1024
+        out, _ = capture_output(
+            engine.simulate,
+            demographic_model=model,
+            contig=contig,
+            samples=samples,
+            slim_script=True,
+            slim_scaling_factor=10,
+            seed=seed,
+        )
+        out = re.sub(
+            'defineConstant\\("trees_file.+;',
+            f'defineConstant("trees_file", "{treefile}");',
+            out,
+        )
+        with open(scriptfile, "w") as f:
+            f.write(out)
+        engine._run_slim(scriptfile, seed=seed)
+        ts_external = tskit.load(treefile)
+        ts_external = engine.recap_and_rescale(
+            ts_external,
+            demographic_model=model,
+            contig=contig,
+            samples=samples,
+            slim_scaling_factor=10,
+            seed=seed,
+        )
+        ts_internal = engine.simulate(
+            demographic_model=model,
+            contig=contig,
+            samples=samples,
+            slim_scaling_factor=10,
+            seed=seed,
+        )
+        tables1 = ts_external.dump_tables()
+        tables2 = ts_internal.dump_tables()
+        assert tables1.nodes == tables2.nodes
+        assert tables1.edges == tables2.edges
+        assert tables1.mutations == tables2.mutations
+
     def test_assert_min_version(self):
         engine = stdpopsim.get_engine("slim")
         with mock.patch(