Skip to content

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

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

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