Minimal updates for msprime 1.0

This is a transitional update, which gets most features working with
msprime 1.0, without starting the process of modernising the catalog
models.

Closes popsim-consortium#748
Closes popsim-consortium#740
Closes popsim-consortium#535

stdpopsim/catalog/BosTau/__init__.py

@@ -115,7 +115,7 @@ def _HolsteinFriesan_1M13():
     inspection of the plots.
     """
     populations = [
-        stdpopsim.Population(id="Holstein-Friesian", description="Holstein-Friesian"),
+        stdpopsim.Population(id="Holstein_Friesian", description="Holstein Friesian"),
     ]
     citations = [_MacLeodEtAl.because(stdpopsim.CiteReason.DEM_MODEL)]
 

stdpopsim/cli.py

@@ -577,17 +577,17 @@ def write_simulation_summary(engine, model, contig, samples, seed=None):
     dry_run_text += f"{indent}Population: number_samples (sampling_time_generations):\n"
     sample_counts = [0] * model.num_sampling_populations
     for s in samples:
-        sample_counts[s.population] += 1
+        sample_counts[s.population] += s.num_samples
     for p in range(0, model.num_sampling_populations):
         pop_name = model.populations[p].id
         sample_time = model.populations[p].sampling_time
         dry_run_text += f"{indent * 2}{pop_name}: "
         dry_run_text += f"{sample_counts[p]} ({sample_time})\n"
     # Get information about relevant contig
     gmap = "None" if contig.genetic_map is None else contig.genetic_map.id
-    mean_recomb_rate = contig.recombination_map.mean_recombination_rate
+    mean_recomb_rate = contig.recombination_map.mean_rate
     mut_rate = contig.mutation_rate
-    contig_len = contig.recombination_map.get_length()
+    contig_len = contig.recombination_map.sequence_length
     dry_run_text += "Contig Description:\n"
     dry_run_text += f"{indent}Contig length: {contig_len}\n"
     dry_run_text += f"{indent}Mean recombination rate: {mean_recomb_rate}\n"

stdpopsim/engines.py

@@ -4,6 +4,7 @@
 import attr
 import msprime
 import stdpopsim
+import numpy as np
 
 logger = logging.getLogger(__name__)
 
@@ -173,16 +174,14 @@ def simulate(
             if msprime_model in self.model_citations:
                 self.citations.extend(self.model_citations[msprime_model])
 
-        demographic_events = demographic_model.demographic_events.copy()
         if msprime_change_model is not None:
+            msprime_model = [msprime_model]
             for t, model in msprime_change_model:
                 if model not in self.supported_models:
                     raise ValueError(f"Unrecognised model '{model}'")
-                model_change = msprime.SimulationModelChange(t, model)
-                demographic_events.append(model_change)
+                msprime_model.append((t, model))
                 if model in self.model_citations:
                     self.citations.extend(self.model_citations[model])
-            demographic_events.sort(key=lambda x: x.time)
 
         if "random_seed" in kwargs.keys():
             if seed is None:
@@ -194,18 +193,25 @@ def simulate(
         # TODO: remove this after a release or two. See #745.
         self._warn_zigzag(demographic_model)
 
-        ts = msprime.simulate(
+        rng = np.random.default_rng(seed)
+        seeds = rng.integers(1, 2 ** 31 - 1, size=2)
+
+        ts = msprime.sim_ancestry(
             samples=samples,
-            recombination_map=contig.recombination_map,
-            mutation_rate=contig.mutation_rate,
-            population_configurations=demographic_model.population_configurations,
-            migration_matrix=demographic_model.migration_matrix,
-            demographic_events=demographic_events,
-            random_seed=seed,
+            recombination_rate=contig.recombination_map,
+            demography=demographic_model.model,
+            ploidy=2,
+            random_seed=seeds[0],
             model=msprime_model,
             end_time=0 if dry_run else None,
             **kwargs,
         )
+        ts = msprime.sim_mutations(
+            ts,
+            end_time=0 if dry_run else None,
+            random_seed=seeds[1],
+            rate=contig.mutation_rate,
+        )
 
         if contig.inclusion_mask is not None:
             ts = stdpopsim.utils.mask_tree_sequence(ts, contig.inclusion_mask, False)

stdpopsim/genetic_maps.py

@@ -4,6 +4,7 @@
 import warnings
 
 import msprime
+import numpy as np
 
 import stdpopsim
 
@@ -88,8 +89,8 @@ def get_chromosome_map(self, id):
         :param str id: The chromosome identifier.
              A complete list of chromosome IDs for each species can be found in the
              "Genome" subsection for the species in the :ref:`sec_catalog`.
-        :rtype: :class:`msprime.RecombinationMap`
-        :return: A :class:`msprime.RecombinationMap` object.
+        :rtype: :class:`msprime.RateMap`
+        :return: A :class:`msprime.RateMap` object.
         """
         chrom = self.species.genome.get_chromosome(id)
         if not self.is_cached():
@@ -100,22 +101,20 @@ def get_chromosome_map(self, id):
         # needs to be redownloaded.
         map_file = self.map_cache_dir / self.file_pattern.format(id=chrom.id)
         if map_file.exists():
-            recomb_map = msprime.RecombinationMap.read_hapmap(str(map_file))
+            recomb_map = msprime.RateMap.read_hapmap(map_file)
         else:
             warnings.warn(
                 "Recombination map not found for chromosome: '{}'"
                 " on map: '{}', substituting a flat map with chromosome "
                 "recombination rate {}".format(id, self.id, chrom.recombination_rate)
             )
-            recomb_map = msprime.RecombinationMap.uniform_map(
-                chrom.length, chrom.recombination_rate
-            )
-        map_length = recomb_map.get_sequence_length()
+            recomb_map = msprime.RateMap.uniform(chrom.length, chrom.recombination_rate)
+        map_length = recomb_map.sequence_length
         if map_length < chrom.length:
             # Extend map to the end of the chromosome.
-            positions = recomb_map.get_positions() + [chrom.length]
-            rates = recomb_map.get_rates() + [0]
-            recomb_map = msprime.RecombinationMap(positions, rates)
+            positions = np.append(recomb_map.position, chrom.length)
+            rates = np.append(recomb_map.rate, 0)
+            recomb_map = msprime.RateMap(positions, rates)
         elif map_length > chrom.length:
             warnings.warn(
                 f"Recombination map has length {map_length}, which is longer than"

stdpopsim/genomes.py

@@ -117,10 +117,8 @@ class Contig:
     :ivar mutation_rate: The rate of mutation per base per generation.
     :vartype mutation_rate: float
     :ivar recombination_map: The recombination map for the region. See the
-        `msprime documentation
-        <https://msprime.readthedocs.io/en/stable/api.html#msprime.RecombinationMap>`_
-        for more details.
-    :vartype recombination_map: msprime.simulations.RecombinationMap
+        :class:`msprime.RateMap` for details.
+    :vartype recombination_map: msprime.RateMap
     :ivar mask_intervals: Intervals to keep in simulated tree sequence, as a list
         of (left_position, right_position), such that intervals are non-overlapping
         and in ascending order. Should have shape Nx2, where N is the number of
@@ -157,8 +155,8 @@ def __str__(self):
             "Contig(length={:.2G}, recombination_rate={:.2G}, "
             "mutation_rate={:.2G}, genetic_map={})"
         ).format(
-            self.recombination_map.get_length(),
-            self.recombination_map.mean_recombination_rate,
+            self.recombination_map.sequence_length,
+            self.recombination_map.mean_rate,
             self.mutation_rate,
             gmap,
         )