PolynomialChaosExpansion: TypeError: only 0-dimensional arrays can be converted to Python scalars

[asher-m]

This was a bug found by Eli a few months ago. As I'm in the process of migrating UncertainSCI to numpy>=2, I thought this would get collected and fixed in the process, but this is turning out to be a bit more complicated than I anticipated.

Here's the MWE Eli provided:

import numpy as np
from UncertainSCI.distributions import BetaDistribution
from UncertainSCI.pce import PolynomialChaosExpansion

beta1 = BetaDistribution(alpha=2.0, beta=5.0, domain=np.array([0.0, 1.0]))
beta2 = BetaDistribution(alpha=3.0, beta=4.0, domain=np.array([0.0, 1.0]))
pce = PolynomialChaosExpansion(
    distribution=[beta1, beta2],
    order=3,
    plabels=["test_param1", "test_param2"],
)

pce.generate_samples()

which yields:

Console Log

Traceback (most recent call last):
  File "/home/elichild01/src/koi-utils/scratch/uncertainscibug.py", line 13, in <module>
    pce.generate_samples()
    ~~~~~~~~~~~~~~~~~~~~^^
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/pce.py", line 222, in generate_samples
    p_standard = self.distribution.polys.wafp_sampling(
                  self.index_set.get_indices(), **self.sampling_options)
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/opolynd.py", line 220, in wafp_sampling
    x = self.idist_mixture_sampling(K, indices, weights=weights,
                                    fast_sampler=fast_sampler)
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/opolynd.py", line 164, in idist_mixture_sampling
    x[:, qd] = idistinv(np.random.random(M), Lambdas[:, qd])
              ~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/families.py", line 463, in fidistinv
    data = self.fidistinv_jacobi_setup(max(n[:]), data)
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/families.py", line 428, in fidistinv_jacobi_setup
    data.append(fidistinv_setup_helper2(ug, idistinv, exponents, 10, self.alpha, self.beta))  # , E_n?
                ~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/families.py", line 189, in fidistinv_setup_helper2
    xgrid[:, q] = idistinv(ugrid[:, q])
                  ~~~~~~~~^^^^^^^^^^^^^
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/families.py", line 426, in idistinv
    return self.idistinv(u, nn)
          ~~~~~~~~~~~~~^^^^^^^
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/families.py", line 364, in idistinv
    x = idistinv_driver(u, n, primitive, ab, supp)
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/opoly1d.py", line 322, in idistinv_driver
    x[j] = optimize.bisect(lambda xx: primitive(xx) - u[j], intervals[j, 0], intervals[j, 1])
          ~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/scipy/optimize/_zeros_py.py", line 595, in bisect
    r = _zeros._bisect(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/scipy/optimize/_zeros_py.py", line 94, in f_raise
    fx = f(x, *args)
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/opoly1d.py", line 322, in <lambda>
    x[j] = optimize.bisect(lambda xx: primitive(xx) - u[j], intervals[j, 0], intervals[j, 1])
                                      ~~~~~~~~~^^^^
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/families.py", line 361, in primitive
    return self.idist(xx, n, M=M)
          ~~~~~~~~~~^^^^^^^^^^^^
  File "/home/elichild01/src/koi-utils/scratch/.venv/lib/python3.14/site-packages/UncertainSCI/families.py", line 335, in idist
    F[np.where(x <= mrs_centroid)] = jacobi_idist_driver(x[np.where(x <= mrs_centroid)], n, self.alpha, self.beta, M)
                                                                    ^^^^^^^^^^^^^^^^^
TypeError: only 0-dimensional arrays can be converted to Python scalars

I'm currently having trouble reproducing this bug on my laptop due to some platform-specific bugs (I think, i.e., I'm getting SystemError: <function> returned a result with exception set; having repro'd this bug on my other machine and searching thus far suggests to me this is not related to UncertainSCI and is instead something with the dependency builds for my laptop, which is itself unrelated but worth investigating). Nonetheless, I verified this on different hardware and can confirm that the bump to numpy>=2 causes this bug to appear.

[asher-m]

The plot thickens: attempting to reproduce this bug on my laptop, I recreated the UncertainSCI python environment with Python 3.9. With Python 3.9 and Numpy 2.0.2 (what pip resolved to satisfy the environment), the bug does not appear.

[asher-m]

Okay, a bit more investigation:

1. The bug does not appear in Numpy 2.0.2 because Numpy 2.0.2 had not yet introduced strict rules about casting to scalar(s). In particular, it is a later version of Numpy that requires explicit casts to scalar to be 0-d arrays.
2. The issue I saw, SystemError: <function> returned a result with an exception set was resulting from the same non-0-d array cast to Python scalar, so there really is no issue with my system or deeper funniness with this bug.

Finally (and this is the real note here): this is going to take a lot more cleanup than I think I'm willing to commit with the test fixes. In particular, this pattern exists in numerous places in families.py, opoly1d.py, and casting.py,

if isinstance(u, float) or isinstance(u, int):
    u = np.asarray([u])
else:
    u = np.asarray(u)

which, as far as I can tell, is used to ensure input(s) are at least 1-d. If my interpretation is correct, then the code should use np.atleast_1d.

However, in many places, this input is then converted to [a] scalar anyways, so it's unclear to me if what we want is a scalar or an array that is at least 1-d. In either case, it seems there is no fundamental need to switch between the two. Choose one and stick with it.

Point is: this cleanup is deeper than I thought and warrants its own focus rather than being included in other things.

I'm going to return my focus to the test fixes for now. I'll make this its own PR when I get to it.

[asher-m]

Okay, one more thing before fully disengaging from this issue: families.py L347 really looks as though functions using this "at least 1-d" pattern expect the input to, in fact, be at most 1-d, e.g. that line and surrounding context reads,

def idistinv(self, u, n):

    if isinstance(u, float) or isinstance(u, int):
        u = np.asarray([u])
    else:
        u = np.asarray(u)
    
    x = np.zeros(u.size)
    (x is assigned later)

In particular, creating x like

x = np.zeros(u.size)

results in an array x that is exactly 1-d.

My suggestion would be to do something like,

def idistinv(self, u, n):
    u = np.atleast_1d(u)
    x = np.empty_like(u)