Merge pull request #549 from grahamgower/Q-warn

Warn when using --slim-scaling-factor with Q!=1.

stdpopsim/__init__.py

@@ -15,6 +15,7 @@
 from . cache import *  # NOQA
 from . citations import *  # NOQA
 from . engines import *  # NOQA
+from . warning_categories import *  # NOQA
 
 # Add imports for all defined species here.
 # We import these here to build the catalog, but the internal functions

stdpopsim/cli.py

@@ -444,7 +444,7 @@ def run_simulation(args):
         write_simulation_summary(engine=engine, model=model, contig=contig,
                                  samples=samples, seed=args.seed)
         if not qc_complete:
-            warnings.warn(
+            warnings.warn(stdpopsim.QCMissingWarning(
                     f"{model.id} has not been QCed. Use at your own risk! "
                     "Demographic models that have not undergone stdpopsim's "
                     "Quality Control procedure may contain implementation "
@@ -453,7 +453,7 @@ def run_simulation(args):
                     "More information about the QC process can be found in "
                     "the developer documentation. "
                     "https://stdpopsim.readthedocs.io/en/latest/development.html"
-                    "#demographic-model-review-process")
+                    "#demographic-model-review-process"))
         ts = engine.simulate(**kwargs)
         summarise_usage()
         if ts is not None:

stdpopsim/slim_engine.py

@@ -596,11 +596,11 @@ def matrix2str(matrix, row_comments=None, col_comment=None, indent=2,
         if count % 2 != 0:
             pop_id = pop_names[pop]
             gen = time / demographic_model.generation_time
-            warnings.warn(
+            warnings.warn(stdpopsim.SLiMOddSampleWarning(
                     f"SLiM simulates diploid individuals, so {n_inds} "
                     f"individuals will be sampled for the {count} haploids "
                     f"requested from population {pop_id} at time {gen}. "
-                    "See #464.")
+                    "See #464."))
         sampling_episodes.append((pop, n_inds, time))
 
     printsc('    // One row for each sampling episode.')
@@ -674,6 +674,14 @@ def simulate(
         if slim_burn_in < 0:
             raise ValueError("slim_burn_in must be non-negative")
 
+        if slim_scaling_factor != 1:
+            warnings.warn(stdpopsim.SLiMScalingFactorWarning(
+                f"You're using a scaling factor ({slim_scaling_factor}). "
+                "This should give similar results for many situations, "
+                "but is not equivalent, especially in the presence of selection. "
+                "When using rescaling, you should be careful---do checks and "
+                "compare results across different values of the scaling factor."))
+
         run_slim = not slim_script
 
         mutation_rate = contig.mutation_rate
@@ -748,7 +756,8 @@ def _run_slim(self, script_file, slim_path=None, seed=None, dry_run=False):
                 if line.startswith("ERROR: "):
                     logger.error(line[len("ERROR: "):])
                 elif line.startswith("WARNING: "):
-                    warnings.warn(line[len("WARNING: "):])
+                    warnings.warn(stdpopsim.UnspecifiedSLiMWarning(
+                        line[len("WARNING: "):]))
                 else:
                     # filter `dbg` function calls that generate output
                     line = line.replace("dbg(self.source); ", "")

stdpopsim/species.py

@@ -8,7 +8,7 @@
 import attr
 import msprime
 
-from . import genomes
+import stdpopsim
 
 logger = logging.getLogger(__name__)
 
@@ -137,10 +137,10 @@ def get_contig(self, chromosome, genetic_map=None, length_multiplier=1):
         # TODO: add non-autosomal support
         if (chromosome is not None and
                 chromosome.lower() in ("x", "y", "m", "mt", "chrx", "chry", "chrm")):
-            warnings.warn(
+            warnings.warn(stdpopsim.NonAutosomalWarning(
                     "Non-autosomal simulations are not yet supported. See "
                     "https://github.com/popsim-consortium/stdpopsim/issues/383 and "
-                    "https://github.com/popsim-consortium/stdpopsim/issues/406")
+                    "https://github.com/popsim-consortium/stdpopsim/issues/406"))
         chrom = self.genome.get_chromosome(chromosome)
         if genetic_map is None:
             logger.debug(f"Making flat chromosome {length_multiplier} * {chrom.id}")
@@ -154,7 +154,7 @@ def get_contig(self, chromosome, genetic_map=None, length_multiplier=1):
             gm = self.get_genetic_map(genetic_map)
             recomb_map = gm.get_chromosome_map(chrom.id)
 
-        ret = genomes.Contig(
+        ret = stdpopsim.Contig(
             recombination_map=recomb_map, mutation_rate=chrom.mutation_rate,
             genetic_map=gm)
         return ret

stdpopsim/warning_categories.py

@@ -0,0 +1,24 @@
+"""
+Specific warning categories to more easily test whether specific warnings
+have been emitted.
+"""
+
+
+class NonAutosomalWarning(UserWarning):
+    pass
+
+
+class QCMissingWarning(UserWarning):
+    pass
+
+
+class SLiMScalingFactorWarning(UserWarning):
+    pass
+
+
+class SLiMOddSampleWarning(UserWarning):
+    pass
+
+
+class UnspecifiedSLiMWarning(UserWarning):
+    pass

tests/test_cli.py

@@ -811,9 +811,8 @@ class TestNonAutosomal(unittest.TestCase):
     # https://github.com/popsim-consortium/stdpopsim/issues/383
     def test_chrX_gives_a_warning(self):
         cmd = "HomSap -D -c chrX -o /dev/null 10".split()
-        with mock.patch("warnings.warn", autospec=True) as mock_warning:
+        with self.assertWarns(stdpopsim.NonAutosomalWarning):
             capture_output(stdpopsim.cli.stdpopsim_main, cmd)
-        mock_warning.assert_called_once()
 
 
 class TestNoQCWarning(unittest.TestCase):
@@ -842,9 +841,8 @@ def tearDown(self):
         self.species.demographic_models.remove(self.model)
 
     def verify_noQC_warning(self, cmd):
-        with mock.patch("warnings.warn", autospec=True) as mock_warning:
+        with self.assertWarns(stdpopsim.QCMissingWarning):
             capture_output(stdpopsim.cli.stdpopsim_main, cmd.split())
-        mock_warning.assert_called_once()
 
     def test_noQC_warning(self):
         self.verify_noQC_warning("EscCol -d FakeModel -D 10")

tests/test_slim_engine.py

@@ -241,21 +241,21 @@ class TestWarningsAndErrors(unittest.TestCase):
     """
     Checks that warning messages are printed when appropriate.
     """
-    def test_odd_sample_warning(self):
+    # this is an expected failure, because no warning should be emitted
+    @unittest.expectedFailure
+    def test_odd_sample_warning_for_even_samples(self):
         cmd = "-q -e slim --slim-script HomSap -d OutOfAfrica_2T12 4 6".split()
-        with mock.patch("warnings.warn", autospec=True) as mock_warning:
+        with self.assertWarns(stdpopsim.SLiMOddSampleWarning):
             capture_output(stdpopsim.cli.stdpopsim_main, cmd)
-        self.assertEqual(mock_warning.call_count, 0)
 
+    def test_odd_sample_warning(self):
         cmd = "-q -e slim --slim-script HomSap -d OutOfAfrica_2T12 4 5".split()
-        with mock.patch("warnings.warn", autospec=True) as mock_warning:
+        with self.assertWarns(stdpopsim.SLiMOddSampleWarning):
             capture_output(stdpopsim.cli.stdpopsim_main, cmd)
-        self.assertEqual(mock_warning.call_count, 1)
 
         cmd = "-q -e slim --slim-script HomSap -d OutOfAfrica_2T12 3 5".split()
-        with mock.patch("warnings.warn", autospec=True) as mock_warning:
+        with self.assertWarns(stdpopsim.SLiMOddSampleWarning):
             capture_output(stdpopsim.cli.stdpopsim_main, cmd)
-        self.assertEqual(mock_warning.call_count, 2)
 
     def triplet(self):
         engine = stdpopsim.get_engine("slim")
@@ -268,11 +268,10 @@ def test_bad_population_size_addSubPop(self):
         model = stdpopsim.PiecewiseConstantSize(100)
         samples = model.get_samples(2)
 
-        with mock.patch("warnings.warn", autospec=True) as mock_warning:
+        with self.assertWarns(stdpopsim.UnspecifiedSLiMWarning):
             engine.simulate(
                     demographic_model=model, contig=contig, samples=samples,
                     slim_scaling_factor=10, dry_run=True)
-        mock_warning.assert_called_once()
 
     def test_no_populations_in_generation_1(self):
         engine, species, contig = self.triplet()
@@ -290,11 +289,10 @@ def test_bad_population_size_addSubpopSplit(self):
                 NA=1000, N1=100, N2=1000, T=1000, M12=0, M21=0)
         samples = model.get_samples(2)
 
-        with mock.patch("warnings.warn", autospec=True) as mock_warning:
+        with self.assertWarns(stdpopsim.UnspecifiedSLiMWarning):
             engine.simulate(
                     demographic_model=model, contig=contig, samples=samples,
                     slim_scaling_factor=10, dry_run=True)
-        mock_warning.assert_called_once()
 
         with self.assertRaises(stdpopsim.SLiMException):
             engine.simulate(
@@ -306,11 +304,10 @@ def test_bad_population_size_setSubpopulationSize(self):
         model = stdpopsim.PiecewiseConstantSize(100, (1000, 1000))
         samples = model.get_samples(2)
 
-        with mock.patch("warnings.warn", autospec=True) as mock_warning:
+        with self.assertWarns(stdpopsim.UnspecifiedSLiMWarning):
             engine.simulate(
                     demographic_model=model, contig=contig, samples=samples,
                     slim_scaling_factor=10, dry_run=True)
-        mock_warning.assert_called_once()
 
         with self.assertRaises(stdpopsim.SLiMException):
             engine.simulate(
@@ -355,11 +352,10 @@ def test_bad_population_size_exp_decline(self):
         model = self.exp_decline()
         samples = model.get_samples(2)
 
-        with mock.patch("warnings.warn", autospec=True) as mock_warning:
+        with self.assertWarns(stdpopsim.UnspecifiedSLiMWarning):
             engine.simulate(
                     demographic_model=model, contig=contig, samples=samples,
                     slim_scaling_factor=10, dry_run=True)
-        mock_warning.assert_called_once()
 
         with self.assertRaises(stdpopsim.SLiMException):
             engine.simulate(
@@ -376,6 +372,26 @@ def test_sample_size_too_big_exp_decline(self):
                     demographic_model=model, contig=contig, samples=samples,
                     slim_scaling_factor=10, dry_run=True)
 
+    # this is an expected failure, because no warning should be emitted
+    @unittest.expectedFailure
+    def test_warning_when_not_scaling(self):
+        with self.assertWarns(stdpopsim.SLiMScalingFactorWarning):
+            for cmd in [
+                    "HomSap 100 -D",
+                    "-e slim HomSap 100 -D",
+                    "-e slim --slim-scaling-factor 1 HomSap 100 -D",
+                    "-e slim --slim-scaling-factor 1.0 HomSap 100 -D",
+                    ]:
+                capture_output(stdpopsim.cli.stdpopsim_main, cmd.split())
+
+    def test_warning_when_scaling(self):
+        for cmd in [
+                "-e slim --slim-scaling-factor 2 HomSap 100 -D",
+                "-e slim --slim-scaling-factor 1000 EscCol 100 -D",
+                ]:
+            with self.assertWarns(stdpopsim.SLiMScalingFactorWarning):
+                capture_output(stdpopsim.cli.stdpopsim_main, cmd.split())
+
 
 class TestSlimAvailable(unittest.TestCase):
     """