MWE for potential SFS issue in stdpopsim

testSFS.py

import stdpopsim
import msprime as msp
import allel
import numpy as np
from stdpopsim import homo_sapiens

chrom = homo_sapiens.genome.chromosomes["chr22"]
recomb_map = chrom.recombination_map()
model = homo_sapiens.TennessenTwoPopOutOfAfrica()
samples = [msp.Sample(population=0,time=0)]*10
# simulate with map
ts = msp.simulate(samples=samples,recombination_map=chrom.recombination_map(),
                  mutation_rate=chrom.default_mutation_rate,**model.asdict())
s = allel.sfs(allel.HaplotypeArray(ts.genotype_matrix()).count_alleles()[:,1])
print("SFS with genetic map")
print(s/np.sum(s))
print("----------------------\n")
# ld prune
haps = allel.HaplotypeArray(ts.genotype_matrix())
ul = allel.locate_unlinked(haps.to_genotypes(ploidy=2).to_n_alt())
s = allel.sfs(haps[ul,:].count_alleles()[:,1])
print("SFS with genetic map, LD pruned")
print(s/np.sum(s))
print("----------------------\n")

# simulate with constant rate
ts = msp.simulate(samples=samples,
                  recombination_rate=chrom.default_recombination_rate,
                  length=chrom.length,
                  mutation_rate=chrom.default_mutation_rate,**model.asdict())
s = allel.sfs(allel.HaplotypeArray(ts.genotype_matrix()).count_alleles()[:,1])
print("SFS with constant rate")
print(s/np.sum(s))
print("----------------------\n")

# ld prune
haps = allel.HaplotypeArray(ts.genotype_matrix())
ul = allel.locate_unlinked(haps)
s = allel.sfs(haps[ul,:].count_alleles()[:,1])
print("SFS with constant rate, LD pruned")
print(s/np.sum(s))