Skip to content

Python bindings silently misread non-C-contiguous matrices (e.g. pandas .values), corrupting the likelihood #1

Description

@colganwi

Summary

The Python bindings silently misread non-C-contiguous NumPy arrays. pylaml.compute_likelihood, pylaml.optimize, and pylaml.topology_search all convert their input matrices in python/src/numpy_conversion.cpp, and those converters assume a C-contiguous (row-major) memory layout:

// numpy_to_character_matrix
auto ptr = static_cast<int32_t*>(buf.ptr);
for (size_t i = 0; i < num_leaves; i++)
    for (size_t j = 0; j < num_chars; j++)
        result[i][j] = ptr[i * num_chars + j];   // hardcoded row-major strides

If the array passed in is not C-contiguous — e.g. a transposed view, a strided slice, or a pandas DataFrame.values (which is frequently Fortran-ordered) — the raw pointer arithmetic reads the wrong elements and the character matrix is effectively transposed/scrambled. No error is raised; the likelihood is just wrong.

Impact

Because the data no longer corresponds to the correct leaves, the likelihood stops rewarding the true tree:

  • the true (generating) tree scores lower than arbitrary/incorrect trees;
  • topology_search "optimizes" a corrupted objective and degrades a good starting tree instead of improving it;
  • randomly permuting the cell→leaf assignment can increase the likelihood (impossible for a correct phylogenetic likelihood);
  • the joint log-likelihood is not equal to the sum of per-character log-likelihoods, and depends on column order.

This is easy to hit in practice: the natural way to feed a character matrix from pandas (df.values) yields an F-contiguous array.

The same bug is present in numpy_to_observation_matrix and numpy_to_mutation_priors.

Reproducer

import numpy as np, pylaml

# Balanced 8-leaf tree
edges = [(8,0),(8,1),(9,2),(9,3),(10,4),(10,5),(11,6),(11,7),
         (12,8),(12,9),(13,10),(13,11),(14,12),(14,13)]
tree = pylaml.make_tree(edges=edges, branch_lengths=[0.4]*14+[0.0], num_leaves=8, root=14)

rng = np.random.RandomState(0)
cm = rng.randint(0, 6, size=(8, 12)).astype(np.int32)

c_order = np.ascontiguousarray(cm)
f_order = np.asfortranarray(cm)          # same logical data, F-contiguous

print(pylaml.compute_likelihood(tree=tree, character_matrix=c_order, nu=0.1, phi=0.05))
print(pylaml.compute_likelihood(tree=tree, character_matrix=f_order, nu=0.1, phi=0.05))
# -> the two values differ; they must be identical

Likewise, for a matrix with heterogeneous per-character alphabet sizes the joint likelihood does not equal the sum over characters and is not invariant to column order.

Root cause

numpy_to_character_matrix, numpy_to_observation_matrix, and numpy_to_mutation_priors in python/src/numpy_conversion.cpp index the raw buffer assuming the default C-contiguous strides instead of honoring buf.strides.

Suggested fix

Read the buffer using its actual strides (or force a C-contiguous copy at the binding boundary). A stride-aware read fixes all entry points without copying.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions