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.

Also updates the DemographicModel class to remove instance variables
that were tied to the old representation and their documentation.

Closes #748
Closes #740
Closes #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,
         )