Python 3.14/NumPy 2.0+ compatibility fixes

dinf/dinf_model.py

@@ -14,9 +14,10 @@
 from .misc import pytree_equal, pytree_shape, pytree_dtype
 
 FeatureCollection = Union[np.ndarray, Dict[str, np.ndarray]]
-FeatureCollection.__doc__ = (
-    """A feature array or a labelled collection of feature arrays."""
-)
+#This bit is commented out to be compatible with Python 3.14:
+#FeatureCollection.__doc__ = (
+#    """A feature array or a labelled collection of feature arrays."""
+#)
 
 
 def _sim_shim(args, *, func, keys):

dinf/feature_extractor.py

@@ -714,8 +714,11 @@ def from_vcf(
             require_phased=self._global_phased,
             rng=rng,
         )
+        sample_indices = {
+            pop: [vb.sample_index_map[id] for id in vb.samples[pop]]
+            for pop in list(vb.samples)
+        }
         num_samples = [len(v) for v in vb.samples.values()]
-        offsets = np.cumsum([0] + num_samples[:-1])
 
         labelled_indexes = {}
         ac0 = np.zeros(len(positions))
@@ -728,8 +731,8 @@ def from_vcf(
                 )
 
             # Subsample individuals.
-            idx = offsets[j] + rng.choice(
-                num_samples[j], size=self._num_individuals[label], replace=False
+            idx = rng.choice(
+                sample_indices[label], size=self._num_individuals[label], replace=False
             )
             labelled_indexes[label] = idx
             H = G[:, idx, : self._ploidy[label]]

dinf/plot.py

@@ -1189,7 +1189,7 @@ def __call__(self, args: argparse.Namespace):
                 ax.set_xlabel("Iteration")
                 ax.xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(integer=True))
 
-                if parameters is not None and x_param != "_Pr":
+                if parameters is not None and x_param != "_Pr" and parameters[x_param].truth is not None:
                     truth = parameters[x_param].truth
                     ax.axhline(truth, c="red", ls="-", alpha=0.7)
 
@@ -1200,7 +1200,7 @@ def __call__(self, args: argparse.Namespace):
 
 def main(args_list=None):
     top_parser = argparse.ArgumentParser(
-        prog="dinf.plot", description="Dinf plotting tools."
+        prog="dinf-plot", description="Dinf plotting tools."
     )
     top_parser.add_argument(
         "-V", "--version", action="version", version=dinf.__version__

dinf/tabulate.py

@@ -217,7 +217,7 @@ def __call__(self, args: argparse.Namespace):
 
 def main(args_list=None):
     top_parser = argparse.ArgumentParser(
-        prog="dinf.tabulate", description="Tabulate Dinf output."
+        prog="dinf-tabulate", description="Tabulate Dinf output."
     )
     top_parser.add_argument(
         "-V", "--version", action="version", version=dinf.__version__

dinf/vcf.py

@@ -24,7 +24,7 @@ def get_samples_from_1kgp_metadata(filename: str, /, *, populations: list) -> di
     :return:
         A dict mapping population names to a list of sample IDs.
     """
-    data = np.recfromtxt(filename, names=True, encoding="ascii")
+    data = np.genfromtxt(filename, names=True, dtype=['<U16', '<U16', '<U16', '<U16', 'i2', '<U16', '<U16'], encoding="ascii").view(np.recarray)
     # Remove related individuals.
     data = data[data.FatherID == "0"]
     data = data[data.MotherID == "0"]
@@ -121,14 +121,24 @@ def get_genotype_matrix(
     """
     gt_list = []
     positions_list = []
+    n_empty = 0
     for variant in vcf(f"{chrom}:{start}-{end}"):
         # XXX: Many variant fields are broken for ploidy > 2.
         # https://github.com/brentp/cyvcf2/issues/227
 
         if not variant.is_snp:
             continue
 
-        a = variant.genotype.array()
+        try:
+            a = variant.genotype.array()
+        except AttributeError as e:
+            n_empty += 1
+            if n_empty >= 10:
+                raise ValueError(
+                    f"At least 10 variants from region {chrom}:{start:d}-{end:d} failed to return valid genotypes, cannot proceed. Did any samples from your metadata get detected in the VCFs?"
+                ) from e
+            else:
+                continue
         gt = a[:, :-1]
 
         # Check phasing. For ploidy == 1, genotypes are reported as unphased.
@@ -209,6 +219,17 @@ class BagOfVcf(collections.abc.Mapping):
     that all individuals in the VCF will be treated as exchangeable.
     """
 
+    sample_index_map: collections.abc.Mapping[str, int]
+    """
+    A dictionary that maps an individual name to a genotype matrix column index.
+
+    The individual names correspond to the VCF columns for which genotypes
+    will be sampled, and the genotype matrix column indices correspond to
+    the genotype matrix containing only the samples requested. This map
+    enables random subsampling of individuals from a population where the
+    individuals may not be in contiguous population blocks in the VCF.
+    """
+
     def __init__(
         self,
         files: Iterable[str | pathlib.Path],
@@ -311,6 +332,10 @@ def _fill_bag(
                         f"Requested individuals not found in {file}: "
                         f"{', '.join(individuals_not_found)}."
                     )
+            #This gets overwritten several times, but that's okay:
+            self.sample_index_map = {
+                id: i for i, id in enumerate(vcf.samples)
+            }
             if contigs is None:
                 try:
                     vcf.seqlens
@@ -360,6 +385,9 @@ def _fill_bag(
 
         # The VCF objects are not valid in child processes, so we set a flag
         # before forking to ensure they're reopened when used in the child.
+        #Note: Given 'spawn' is used instead of 'fork', this doesn't actually
+        # have any impact, but could be part of the code path if switched
+        # back to 'fork'.
         self._needs_reopen = False
         os.register_at_fork(before=self._close)
 

examples/gutenkunst2009/model.py

@@ -14,7 +14,7 @@
 )
 contig_lengths = dinf.get_contig_lengths(
     "GRCh38_full_analysis_set_plus_decoy_hla.fa.fai",
-    keep_contigs={f"chr{c + 1}" for c in range(21)},  # Exclude chrX, etc.
+    keep_contigs={f"chr{c + 1}" for c in range(22)},  # Exclude chrX, etc.
 )
 num_individuals = 64
 recombination_rate = 1.25e-8

examples/one_population_growth/model.py

@@ -12,11 +12,11 @@
 
 population_name = "CEU"
 samples = dinf.get_samples_from_1kgp_metadata(
-    "20130606_g1k_3202_samples_ped_population.txt", [population_name]
+    "20130606_g1k_3202_samples_ped_population.txt", populations=[population_name]
 )
 contig_lengths = dinf.get_contig_lengths(
     "GRCh38_full_analysis_set_plus_decoy_hla.fa.fai",
-    keep_contigs={f"chr{c + 1}" for c in range(21)},  # Exclude chrX, etc.
+    keep_contigs={f"chr{c + 1}" for c in range(22)},  # Exclude chrX, etc.
 )
 recombination_rate = 1.25e-8
 mutation_rate = 1.25e-8

examples/post_ooa/model.py

@@ -12,11 +12,11 @@
 
 populations = ["CEU", "CHB"]
 samples = dinf.get_samples_from_1kgp_metadata(
-    "20130606_g1k_3202_samples_ped_population.txt", populations
+    "20130606_g1k_3202_samples_ped_population.txt", populations=populations
 )
 contig_lengths = dinf.get_contig_lengths(
     "GRCh38_full_analysis_set_plus_decoy_hla.fa.fai",
-    keep_contigs={f"chr{c + 1}" for c in range(21)},  # Exclude chrX, etc.
+    keep_contigs={f"chr{c + 1}" for c in range(22)},  # Exclude chrX, etc.
 )
 recombination_rate = 1.25e-8
 mutation_rate = 1.25e-8

examples/two_pop_ooa/model.py

@@ -12,11 +12,11 @@
 
 populations = ["YRI", "CEU"]
 samples = dinf.get_samples_from_1kgp_metadata(
-    "20130606_g1k_3202_samples_ped_population.txt", populations
+    "20130606_g1k_3202_samples_ped_population.txt", populations=populations
 )
 contig_lengths = dinf.get_contig_lengths(
     "GRCh38_full_analysis_set_plus_decoy_hla.fa.fai",
-    keep_contigs={f"chr{c + 1}" for c in range(21)},  # Exclude chrX, etc.
+    keep_contigs={f"chr{c + 1}" for c in range(22)},  # Exclude chrX, etc.
 )
 recombination_rate = 1.25e-8
 mutation_rate = 1.25e-8

tests/conftest.py

@@ -5,7 +5,7 @@
 
 
 @pytest.fixture(scope="session")
-@pytest.mark.usefixtures("tmp_path_factory")
+#@pytest.mark.usefixtures("tmp_path_factory")
 def discriminator_file(tmp_path_factory):
     discriminator_file = tmp_path_factory.mktemp("discr") / "discriminator.nn"
     ex = "examples/bottleneck/model.py"
@@ -28,8 +28,8 @@ def discriminator_file(tmp_path_factory):
 
 
 @pytest.fixture(scope="session")
-@pytest.mark.usefixtures("tmp_path_factory")
-@pytest.mark.usefixtures("discriminator_file")
+#@pytest.mark.usefixtures("tmp_path_factory")
+#@pytest.mark.usefixtures("discriminator_file")
 def data_file(tmp_path_factory, discriminator_file):
     data_file = tmp_path_factory.mktemp("data") / "data.npz"
     ex = "examples/bottleneck/model.py"
@@ -51,7 +51,7 @@ def data_file(tmp_path_factory, discriminator_file):
 
 
 @pytest.fixture(scope="session")
-@pytest.mark.usefixtures("tmp_path_factory")
+#@pytest.mark.usefixtures("tmp_path_factory")
 def mc_outdir(tmp_path_factory):
     output_folder = tmp_path_factory.mktemp("out")
     ex = "examples/bottleneck/model.py"

tests/test_vcf.py

@@ -19,6 +19,10 @@
 FamilyID SampleID FatherID MotherID Sex Population Superpopulation
 HG00096 HG00096 0 0 1 GBR EUR
 HG00097 HG00097 0 0 2 GBR EUR
+HG00180 HG00180 0 0 2 FIN EUR
+HG00181 HG00181 0 0 1 FIN EUR
+HG00234 HG00234 0 0 1 GBR EUR
+HG00235 HG00235 0 0 2 GBR EUR
 SH001 HG00403 0 0 1 CHS EAS
 SH001 HG00404 0 0 2 CHS EAS
 SH001 HG00405 HG00403 HG00404 2 CHS EAS
@@ -32,13 +36,15 @@ def test_simple(self, tmp_path):
         with open(filename, "w") as f:
             f.write(_1kgp_test_metadata)
         gbr = dinf.get_samples_from_1kgp_metadata(filename, populations=["GBR"])
-        assert gbr == {"GBR": ["HG00096", "HG00097"]}
+        assert gbr == {"GBR": ["HG00096", "HG00097", "HG00234", "HG00235"]}
+        fin = dinf.get_samples_from_1kgp_metadata(filename, populations=["FIN"])
+        assert fin == {"FIN": ["HG00180", "HG00181"]}
         chs = dinf.get_samples_from_1kgp_metadata(filename, populations=["CHS"])
         assert chs == {"CHS": ["HG00403", "HG00404"]}
         samples = dinf.get_samples_from_1kgp_metadata(
-            filename, populations=["GBR", "CHS"]
+            filename, populations=["GBR", "FIN", "CHS"]
         )
-        assert samples == dict(**gbr, **chs)
+        assert samples == dict(**gbr, **fin, **chs)
 
 
 class TestGetContigLengths: