Skip to content

genetic_value() writes past the node array when the ancestral allele is causal #191

Closed

Description

@gregorgorjanc

With help from an LLM, I found that genetic_value() writes past the node array when the ancestral allele is causal - in that case tree.virtual_root is used, but the output array has space only for real nodes.

Normal execution appears correct because Numba disables bounds checking by default. The existing test passes normally:

pytest -q tests/test_genetic_value.py::TestGenotype::test_allele_freq_one

giving

.                                                     [100%]
1 passed in 8.87s

but exposes the invalid write when bounds checking is enabled:

NUMBA_BOUNDSCHECK=1 pytest -q tests/test_genetic_value.py::TestGenotype::test_allele_freq_one

giving

F                                                     [100%]
========================= FAILURES ==========================
_____________ TestGenotype.test_allele_freq_one _____________

self = <tests.test_genetic_value.TestGenotype object at 0x14a900ef0>

    def test_allele_freq_one(self):
        ts = binary_tree()
        tables = ts.dump_tables()
        tables.sites.add_row(4, "A")
        tables.mutations.add_row(site=4, node=0, derived_state="T")
        tables.mutations.add_row(site=4, node=0, derived_state="A", parent=9)
        ts = tables.tree_sequence()
        trait_df = pd.DataFrame(
            {
                "site_id": [4],
                "effect_size": [1],
                "trait_id": [0],
                "causal_allele": ["A"],
            }
        )
>       genetic_result = tstrait.genetic_value(ts=ts, trait_df=trait_df)
                         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

tests/test_genetic_value.py:454: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
tstrait/genetic_value.py:219: in genetic_value
    genetic_result = genetic._run(level)
                     ^^^^^^^^^^^^^^^^^^^
tstrait/genetic_value.py:133: in _run
    genetic_value = self._node_genetic_values(
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <tstrait.genetic_value._GeneticValue object at 0x14a795d30>
tree = <tskit.trees.Tree object at 0x14a6cb4d0>
site = Site(id=4, position=4.0, ancestral_state='A', mutations=[Mutation(id=9, site=4, node=0, derived_state='T', parent=-1, ...d=10, site=4, node=0, derived_state='A', parent=9, metadata=b'', time=nan, edge=0, inherited_state='T')], metadata=b'')
causal_allele = 'A', effect_size = 1

    def _node_genetic_values(self, tree, site, causal_allele, effect_size):
        """
        Returns a numpy array with node genetic values.
        """
        has_mutation = np.zeros(self.ts.num_nodes + 1, dtype=bool)
        state_transitions = {tree.virtual_root: site.ancestral_state}
        for m in site.mutations:
            state_transitions[m.node] = m.derived_state
            has_mutation[m.node] = True
        stack = numba.typed.List()
        for node, allele in state_transitions.items():
            if allele == causal_allele:
                stack.append(node)
    
        if len(stack) == 0:
            genetic_value = np.zeros(self.ts.num_nodes)
        else:
>           genetic_value = _compute_nodes_genetic_value(
                left_child_array=tree.left_child_array,
                right_sib_array=tree.right_sib_array,
                stack=stack,
                has_mutation=has_mutation,
                num_nodes=self.ts.num_nodes,
                effect_size=effect_size,
            )
E           IndexError: index is out of bounds

tstrait/genetic_value.py:97: IndexError
================== short test summary info ==================
FAILED tests/test_genetic_value.py::TestGenotype::test_allele_freq_one - IndexError: index is out of bounds
1 failed in 2.53s

Not sure how serious this is! The computations seem correct. I don't know where Python writes when it is out of bounds - to another random place in memory?, if so, would be good to fix this.

Activity

  1. gregorgorjanc commented on Sep 1, 2026

    @gregorgorjanc
    MemberAuthor

    Here is a small script

    import pandas as pd
    import tskit
    import tstrait
    
    
    ts = tskit.Tree.generate_balanced(2, span=1).tree_sequence
    tables = ts.dump_tables()
    
    # Assign each sample node to one haploid individual.
    tables.individuals.add_row()
    tables.individuals.add_row()
    tables.nodes.individual = [0, 1, tskit.NULL]
    
    site_id = tables.sites.add_row(position=0, ancestral_state="T")
    tables.mutations.add_row(site=site_id, node=0, derived_state="A")
    ts = tables.tree_sequence()
    
    trait_df = pd.DataFrame(
        {
            "site_id": [site_id],
            "effect_size": [1.0],
            "trait_id": [0],
            "causal_allele": ["T"],
        }
    )
    
    print("num_nodes:", ts.num_nodes)
    print("virtual_root:", ts.first().virtual_root)
    print("tstrait:", tstrait.__file__)
    for level in ["node", "edge", "individual"]:
        print(f"\n{level} values")
        print(tstrait.genetic_value(ts, trait_df, level=level))

    giving

    ❯ mamba run -n arg_work python reproduce.py
    
    num_nodes: 3
    virtual_root: 3
    tstrait: /Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/__init__.py
    
    node values
       trait_id  node_id  genetic_value
    0         0        0            0.0
    1         0        1            1.0
    2         0        2            1.0
    
    edge values
       trait_id  edge_id  genetic_value
    0         0        0            0.0
    1         0        1            1.0
    
    individual values
       trait_id  individual_id  genetic_value
    0         0              0            0.0
    1         0              1            1.0
    
    ❯ mamba run -n arg_work env NUMBA_BOUNDSCHECK=1 python reproduce.py 
    
    num_nodes: 3
    virtual_root: 3
    tstrait: /Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/__init__.py
    
    node values
    
    Traceback (most recent call last):
      File "/Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/reproduce.py", line 32, in <module>
        print(tstrait.genetic_value(ts, trait_df, level=level))
              ~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^
      File "/Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/genetic_value.py", line 219, in genetic_value
        genetic_result = genetic._run(level)
      File "/Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/genetic_value.py", line 133, in _run
        genetic_value = self._node_genetic_values(
            tree=tree,
        ...<2 lines>...
            effect_size=data.effect_size,
        )
      File "/Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/genetic_value.py", line 97, in _node_genetic_values
        genetic_value = _compute_nodes_genetic_value(
            left_child_array=tree.left_child_array,
        ...<4 lines>...
            effect_size=effect_size,
        )
    IndexError: index is out of bounds
    
    ERROR conda.cli.main_run:execute(125): `conda run env NUMBA_BOUNDSCHECK=1 python reproduce.py` failed. (See above for error)
    ❯ mamba run -n arg_work env NUMBA_DISABLE_JIT=1 python reproduce.py
    
    num_nodes: 3
    virtual_root: 3
    tstrait: /Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/__init__.py
    
    node values
    
    Traceback (most recent call last):
      File "/Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/reproduce.py", line 32, in <module>
        print(tstrait.genetic_value(ts, trait_df, level=level))
              ~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^
      File "/Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/genetic_value.py", line 219, in genetic_value
        genetic_result = genetic._run(level)
      File "/Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/genetic_value.py", line 133, in _run
        genetic_value = self._node_genetic_values(
            tree=tree,
        ...<2 lines>...
            effect_size=data.effect_size,
        )
      File "/Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/genetic_value.py", line 97, in _node_genetic_values
        genetic_value = _compute_nodes_genetic_value(
            left_child_array=tree.left_child_array,
        ...<4 lines>...
            effect_size=effect_size,
        )
      File "/Users/ggorjanc/Storages/GitBox/tskit/tstrait/tstrait/tstrait/genetic_value.py", line 25, in _compute_nodes_genetic_value
        genetic_value[parent_node_id] = effect_size
        ~~~~~~~~~~~~~^^^^^^^^^^^^^^^^
    IndexError: index 3 is out of bounds for axis 0 with size 3
    
    ERROR conda.cli.main_run:execute(125): `conda run env NUMBA_DISABLE_JIT=1 python reproduce.py` failed. (See above for error)
    
  2. jeromekelleher commented on Sep 2, 2026

    @jeromekelleher
    Member

    I think this is probably reasonably serious and we should fix it.

  3. jeromekelleher commented on Sep 2, 2026

    @jeromekelleher
    Member

    Great catch by the way!

  4. jeromekelleher commented on Sep 2, 2026

    @jeromekelleher
    Member

    I'll take a look at this now.

  5. jeromekelleher commented on Sep 2, 2026

    @jeromekelleher
    Member

    Closed in #192

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions