Error in sim_mutations: "AttributeError: 'generator' object has no attribute 'tables'"

[alanrogers]

I'm trying to move to the 1.x API for msprime and am getting the following error from sim_mutations:

Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/msprime/mutations.py", line 1325, in sim_mutations
    tables = tree_sequence.tables
AttributeError: 'generator' object has no attribute 'tables'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/rogers/Dropbox/tch/archgen/msp/mwe.py", line 108, in <module>
    sim = msp.sim_mutations(treeseq, rate = u_per_site)
  File "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/msprime/mutations.py", line 1327, in sim_mutations
    raise ValueError("First argument must be a TreeSequence instance.")
ValueError: First argument must be a TreeSequence instance.`

Would be grateful for suggestions. My code is below:

#!/usr/bin/python3

# Before running this: "python3 pip install msprime"
import msprime as msp
import os, sys

# time parameters in generations
Txyn = 25920
Tnd = 15000
Txy = 3788
Td = 1734      # age of Denisova fossil
Ta = 1760      # age of Altai fossil
Talpha = 1897     # time of Neanderthal admixture

# population sizes
Nxyn = 64964.1/2.0  # ancestral population
Nxy = 44869.2/2.0
Nn = 9756.8/2.0
Nx = 20000/2.0 # modern Africa
Ny = 20000/2.0 # modern Europe

# admixture
mAlpha = 0.05

nchromosomes = 1      # number of chromosomes
basepairs = 2e6   # number of nucleotides per chromosome
u_per_site = 1.4e-8 # mutation

c = 1e-8 # recombination rate per base pair per generation

# One haploid sample from each of 3 populations: two modern (X,Y),
# and one archaic (N).
samples = [
    msp.SampleSet(num_samples=1, population=0, time=0, ploidy=1),  # pop X
    msp.SampleSet(num_samples=1, population=1, time=0, ploidy=1),  # pop Y
    msp.SampleSet(num_samples=1, population=2, time=Ta, ploidy=1), # pop N
]

# Population configurations. 
dem = msp.Demography()
dem.add_population(
    name = "X",
    description = "African",
    initial_size = Nx
)
dem.add_population(
    name = "Y",
    description = "Eurasian",
    initial_size = Ny
)
dem.add_population(
    name = "Y2",
    description = "Y before alpha gene flow",
    initial_size = Ny
)
dem.add_population(
    name = "N",
    description = "Neanderthal",
    initial_size = Nn
)
dem.add_population(
    name = "XY",
    description = "early modern humans",
    initial_size = Nxy
)
dem.add_population(
    name = "XYN",
    description = "ancestral population",
    initial_size = Nxyn
)

# N->Y gene flow
dem.add_admixture(
    time = Talpha,
    derived = "Y",
    ancestral = ["Y2", "N"],
    proportions = [1-mAlpha, mAlpha]
)

# X-Y split
dem.add_population_split(
    time = Txy,
    derived = ["X", "Y2"],
    ancestral = "XY"
)

# XY-N split
dem.add_population_split(
    time = Txyn,
    derived = ["XY", "N"],
    ancestral = "XYN"
)
    
# Seed for random number generator. Read from /dev/urandom.
seed = int.from_bytes(os.urandom(4), "big")

# Simulate gene genealogy; return iterator over TreeSequence replicates
treeseq = msp.sim_ancestry(samples = samples,
                           demography = dem,
                           sequence_length = basepairs,
                           random_seed = seed,
                           recombination_rate = c,
                           num_replicates = nchromosomes
)
    
# Put mutations onto the gene genealogy
sim = msp.sim_mutations(treeseq, rate = u_per_site)

[alanrogers]

I think I found an anwer to my own question. When sim_ancestry is called without num_replicates, it returns a value of type tskit.trees.TreeSequence:

>>> ts = msp.sim_ancestry(2, sequence_length=100, random_seed=1234)
>>> type(ts)
<class 'tskit.trees.TreeSequence'>

But when it's called with num_replicates, it returns a value of type generator:

>>> ts = msp.sim_ancestry(2, sequence_length=100, random_seed=1234, num_replicates=2)
>>> type(ts)
<class 'generator'>

sim_mutations wants a value of the former type, not a generator.

[alanrogers]

Responding to my own previous post: It might be better just to remove the num_replicates argument from sim_ancestry. The goal of that argument is to spare the user from writing a "for" loop. But if one still needs to loop across replicates when processing the output, then what has been gained? Why not put sim_ancestry and sim_mutations inside the loop? The API would be simpler, and I think the user's code would be clearer.

[jeromekelleher]

"Performance" is the short answer @alanrogers - there's a substantial perf difference when you're doing lots of small simulations. I think that was the motivation in the first place - I agree it's not a great pattern, but the perf difference is compelling.

[alanrogers]

I would be happy to work on this in collaboration with @mtisla. Should we just clone the repo, create a branch, and start coding, or are there other steps to take first?

[benjeffery]

Our development guide at https://tskit.dev/msprime/docs/stable/development.html has a getting started section. It also links out to the tskit docs where a suggested github workflow is detailed.

[jeromekelleher]

Great, thanks @alanrogers - please do reach out if anything is unclear!