Solve numpy/scipy dependency issues

[jekoorn]

In the end #383 (comment) comes to bite us in the .... .
A newer version of numpy introduces a matching issue between the scipy integration, and the numpy output.

An example error log is the following, when running with the FFNS scheme:

 • F2_charm at 1 pts

                                                                                                                       Calculation
yadism took off! please stay tuned ...
a heavy vv gluon is made, with variation -1
Traceback (most recent call last):
  File "/git/yad-4ihou/run.py", line 150, in <module>
    F2_tot, FL, F2_c, FL_c = runner(
                             ~~~~~~^
        Ebeam=1000,
        ^^^^^^^^^^^
    ...<6 lines>...
        theory_path=f"/git/yad-4ihou/{theory}.yaml"
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/git/yad-4ihou/run.py", line 135, in runner
    F2_c = compute_observable("F2_charm", xgrid, Q2grid, Ebeam, PTO, channel, particle, target,runningscales, PDF, theory_path)
  File "/git/yad-4ihou/run.py", line 77, in compute_observable
    output = yadism.run_yadism(theory, observable)
  File "/git/yad-4ihou/src/yadism/__init__.py", line 31, in run_yadism
    return runner.get_result()
           ~~~~~~~~~~~~~~~~~^^
  File "/git/yad-4ihou/src/yadism/runner.py", line 292, in get_result
    results[idx] = elem.get_result()
                   ~~~~~~~~~~~~~~~^^
  File "/git/yad-4ihou/src/yadism/esf/esf.py", line 189, in get_result
    self.compute_local()
    ~~~~~~~~~~~~~~~~~~^^
  File "/git/yad-4ihou/src/yadism/esf/esf.py", line 143, in compute_local
    val, err = conv.convolve_vector(
               ~~~~~~~~~~~~~~~~~~~~^
        rsl, self.info.configs.managers["interpolator"], convolution_point
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/git/yad-4ihou/src/yadism/esf/conv.py", line 230, in convolve_vector
    c, e = convolution(cf, convolution_point, polynomial_f)
           ~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/git/yad-4ihou/src/yadism/esf/conv.py", line 151, in convolution
    res, err = scipy.integrate.quad(
               ~~~~~~~~~~~~~~~~~~~~^
        quad_ker,
        ^^^^^^^^^
    ...<4 lines>...
        points=breakpoints,
        ^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/miniconda/envs/yadism-d/lib/python3.13/site-packages/scipy/integrate/_quadpack_py.py", line 479, in quad
    retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
                   points)
  File "//miniconda/envs/yadism-d/lib/python3.13/site-packages/scipy/integrate/_quadpack_py.py", line 639, in _quad
    return _quadpack._qagpe(func, a, b, the_points, args, full_output,
           ~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                            epsabs, epsrel, limit)
                            ^^^^^^^^^^^^^^^^^^^^^^
TypeError: only 0-dimensional arrays can be converted to Python scalars

For now in this PR I have locked the numpy version to 2.2.4, but the purpose of this PR is to fix the matching for good.

Root Cause

In src/yadism/coefficient_functions/heavy/n3lo/init.py, interpolator() returns a raw RectBivariateSpline and RectBivariateSpline.__call__ with grid=True (default if unspecified) always returns a 2D (1,1) array, never a bare scalar. Older numpy silently coerced that to a float but numpy >= 2.4 turned that into a hard TypeError.

[jekoorn]

The issue was that coeff_iterpol in the N3LO interpolator sometimes returns a scalar, sometimes an array. Scalarizing the output seems to do the trick. I explicitly checked it under numpy 2.4.6

[felixhekhorn]

sometimes

"sometimes" makes me nervous 🙈 is there no clear reason? e.g. based on the argument types?

you still have poetry.lock in the changed files - are we sure we want that? I'm not sure about the answer myself ...

[jekoorn]

Sometimes is indeed a scary word, but I will have to look into what the interpolator does exactly and what’s inside the N3LO grids before verifying this🙈

we do really have to think about putting in soms proper benchmarks… And the test that I modified in the last dependency pr

[Radonirinaunimi]

We should implement the fix at the core so that all calls (not only f2_nc and fl_nc here) get the correct behavior. See the following patch:

diff --git a/src/yadism/coefficient_functions/heavy/n3lo/__init__.py b/src/yadism/coefficient_functions/heavy/n3lo/__init__.py
index 7f095724..9dbe1180 100644
--- a/src/yadism/coefficient_functions/heavy/n3lo/__init__.py
+++ b/src/yadism/coefficient_functions/heavy/n3lo/__init__.py
@@ -22,7 +22,10 @@ def interpolator(coeff, nf, variation):
     coeff = np.load(grid_path / grid_name)
     grid_interpolator = RectBivariateSpline(xi_grid, eta_grid, coeff)

-    # store result
-    interpolators[grid_name] = grid_interpolator
+    def evaluate(xi, eta):
+        # `grid=False`: point evaluation. `.item()` avoids array->scalar `TypeError`.
+        return grid_interpolator(xi, eta, grid=False).item()

-    return grid_interpolator
+    interpolators[grid_name] = evaluate
+
+    return evaluate

With this, we do not also need to lock numpy.

PS: I updated the description of the PR so that the root cause is explained for future reference.

[jekoorn]

Amazing, thank you @Radonirinaunimi