CovertLab / wcEcoli

Whole Cell Model of E. coli
Other
18 stars 4 forks source link

Enforce integer sampling from multinomial #1399

Open thalassemia opened 1 year ago

thalassemia commented 1 year ago

When numpy.random.RandomState.multinomial is called in an environment with Cython<3.0.0 (currently, we use 0.29.34), it accepts non-integer n values with no problem (source code for reference). However, something changes in Cython>=3.0.0 which causes the multinomial function to raise an error when n is not an integer. In our model, this happens in at least four places: here, here, here, and here. The solution on our end is probably just to cast each of our n arguments to int once we upgrade Cython.

Digging deeper into the issue, I noticed that the TypeError for floating point values is only raised when a Cython function has a parameter of the type np.npy_intp. If the parameter has any other integer type (e.g. int, long, np.npy_int), it is simply cast to int with no errors raised.

import numpy as np
cimport numpy as np

def func(int length):
    # func(1.1) returns 1
    return length

def weird_func(np.npy_intp length):
    # weird_func(1.1) raises TypeError: 'float' object cannot be interpreted as an integer
    return length

Does anyone know if this issue is platform-specific (I tested on Ubuntu 22.04 x86) or intended behavior?

1fish2 commented 1 year ago
thalassemia commented 1 year ago

Good catch. I only see the issue when calling from Python. I have not tried calling from Cython. In my testing, the TypeError appears starting with Numpy 1.26. Since this is the first Numpy release with official support for Cython 3.0.0, I would not be surprised if the Numpy maintainers also used it to build the wheels starting with that release.

1fish2 commented 1 year ago

That explains it. I'm trying to properly install Numpy 1.26.1, which in turn requires building a fresh pyenv with updated Numba and llvmlite and maybe Aesara and scipy??? These things have inadequate release notes and installation docs. (How to use Apple's Accelerate library now? Install with --no-binary numpy? With the right ~/.numpy-site.cfg? Update to macOS Sonoma?).

Numpy updates often require client code changes (e.g. np.product() -> np.prod()), adapting to subtle API semantic changes (1.25.0 no longer treats size=1 arrays as scalars), API deprecations & removals (e.g. fix(y=...), besides new features (matrix multiply a @ b and a @= b).

Many of these changes can't be found by static analysis, so we should at least be alert to fix deprecation warnings. (The Jenkins builds now show warnings.)

1fish2 commented 1 year ago

I think the cause of this API change was a changed definition of np.npy_intp:

  1. Numpy 1.26.0 took over managing their own, version-specific API bindings rather than having them embedded in Cython and subject to version skew.
  2. That shifted the definition of np.npy_intp from int in Cython's UFuncs.pyx to Py_ssize_t in NumPy's __init__.cython-30.pxd.
1fish2 commented 1 year ago

The numpy bug is tagged for the 1.26.2 release milestone. If they fix it by the time we're ready for a new NumPy release, we won't need the int() workarounds although we'll still need the np.product() -> np.prod() update.

thalassemia commented 1 year ago

Thanks for linking the bug report you made! I find this kind of low-level stuff pretty fascinating. I cloned the Numpy repo and edited __init__.cython-30.pxd to define npy_intp as int instead of Py_ssize_t (this was on the main branch for Numpy 2.0). Building and installing allows me to do np.random.multinomial(20.0, [1/6.]*6, size=1) with no error. Changing the line back to Py_ssize_t and reinstalling causes the function call to return a TypeError again.

thalassemia commented 1 year ago

I think I've teased out all the details here. As noted by the Cython maintainer, Py_ssize_t uses an index-specific conversion protocol that is different from that used for other integer types. This protocol calls the PyNumber_Index function here. This PyNumber_Index function calls PyIndex_Check here, and PyIndex_Check enforces the integer type here (tp_as_number->nb_index is a null pointer for float objects).

Interestingly, commenting out the PyNumber_Index line in the generated Cython code and manually recompiling still leads to a TypeError being raised in the subsequent call to PyInt_AsSsize_t. For Python >3.0, PyInt_AsSsize_t is defined to be PyLong_AsSsize_t, which enforces the integer type here.

Conversely, defining a parameter as int generates Cython code that calls __Pyx_PyNumber_IntOrLong here. This function converts float objects to int here (the nb_int slot of float objects is defined here).

This all aligns with what you said about Python only enforcing integer types for indices, exactly what the Py_ssize_t type was designed for. As for why Numpy is using an index type for multinomial, I found a very old bug fix that changed many variables in the mtrand module to npy_intp since they were being used as indices. In the multinomial function, this variable declaration was changed to npy_intp, forcing the input type to change as well to perform this assignment. This is unnecessary because this variable is never used as an index.

Looking at the current Numpy codebase, I can see that this issue only affects the legacy RandomState class and can be fixed by simply updating the method signature here (the variable that n is assigned to has already been reverted to a long here).

1fish2 commented 11 months ago

NumPy 1.26.2 has @thalassemia's bug fix and the broader fix for Py_intptr_t integer conversion.

I'll test it soon.