Fix inconsistent out-of-range evaluation between Tabulated1D scalar and array paths

[kawacukennedy]

Description

Evaluating openmc.data.Tabulated1D on values outside its tabulated range currently gives different answers depending on whether the input is a scalar or an array:

>>> tab = openmc.data.Tabulated1D([1, 2, 3], [4, 5, 6], [], [2])
>>> tab(0)
4
>>> tab([0])
array([0])

The scalar path (_interpolate_scalar) returns the value at the nearest tabulated endpoint, while the array path in __call__ initializes the output with zeros and only fills points that fall inside an interpolation region, leaving out-of-range entries at zero.

This PR makes the array path assign the boundary values (y[0] / y[-1]) to out-of-range points so that both paths agree. This is also consistent with the existing precision handling at the domain edges (np.isclose checks), which already assigns endpoint values to near-edge points.

Changes

• openmc/data/function.py:
  • Tabulated1D.__call__: out-of-range points now receive the value of the nearest tabulated endpoint, matching the scalar path.
  • Class docstring: documented the out-of-range behavior.
  • sum_functions: each tabulated component is now explicitly evaluated only where it is defined (points outside a component's own tabulated range contribute zero). This preserves the existing behavior of combined functions — e.g., FissionEnergyRelease.recoverable, total, and the q_* properties, which combine components that may cover different incident energy ranges on a union grid — independently of the new out-of-range semantics.

Testing

Added tests/unit_tests/test_function.py covering:

• scalar/array agreement across, below, and above the tabulated range (the issue's regression case),
• all five ENDF interpolation schemes,
• multi-region functions,
• exact-endpoint and floating-point-precision edge cases,
• multidimensional input shape preservation,
• sum_functions behavior for components with differing domains and for polynomial+tabulated combinations.

Local results: all 11 new tests pass; existing unit tests that exercise these code paths were compared before/after the change with identical outcomes (failures observed locally are due to no nuclear data being configured and are present on unmodified develop as well).

Fixes: #4041

Signed-off-by: Engineer kawacukent@gmail.com

[CAOShurong]

Context

Reviewed exact head b39ee90 against develop base 86ceaad.

Summary

The scalar/array boundary fix is useful and the new tests cover its intended interpolation behavior. I found one blocking dtype regression in sum_functions(), however: the new accumulator inherits an integer union-grid dtype and cannot add the floating-point values returned by supported functions.

Detailed findings

Blocking issue

• [Major] Please use an accumulator dtype compatible with the evaluated function values and add an integer-grid Tabulated1D + Polynomial regression. On the base, the supported combination returns [10.5, 20.0, 29.5] as float64; this head raises UFuncOutputCastingError while adding float64 into int64.

Verified areas

• Purpose/scope: focused fix for #4041, with no new dependency or public API.
• Correctness/testing: the 11 new tests pass locally; compileall and diff checking pass. The public rollup currently reports all 18 contexts successful.
• Physics/design/performance/docs: no new physics model or transport-loop allocation; the endpoint semantics are documented and the overall design remains localized.

I used AI assistance to help inspect the repository and run the base/head verification; I checked the exact diff, reproducer, and results before submitting this review.

[CAOShurong]

np.zeros_like(x) inherits x's dtype. When a tabulated grid is integer-valued, y is int64, so adding a Polynomial (or any floating-point function result) raises UFuncOutputCastingError. This worked on develop and is used by supported sum_functions() call paths. Please initialize an accumulator that can safely represent the evaluated values and cover an integer-grid Tabulated1D + Polynomial case.

[kawacukennedy]

Thanks @CAOShurong for the careful review and for catching the dtype regression — that is exactly right.

Reproduced: with an integer-valued grid, np.zeros_like(x) inherits int64, so adding a Polynomial result raised UFuncTypeError: Cannot cast ufunc 'add' output from dtype('float64') to dtype('int64') with casting rule 'same_kind'.

Fixed (commit 56a8e5a): the accumulator in sum_functions is now np.zeros_like(x, dtype=float), restoring the float-promoting behavior of the previous sum(f(x) for f in funcs).

Added regression test: test_sum_functions_integer_grid combines Tabulated1D([2, 4], [10, 20]) (integer grid) with Polynomial((1.0, -0.5)) and asserts the result is float64 with the expected values [10.0, 19.0] — the exact case you flagged.

The full tests/unit_tests/test_function.py suite (12 tests) passes locally. Environmental openmc.lib/nuclear-data errors observed in neighboring test files are pre-existing on unmodified develop and unrelated to these Python-only changes.

[dylanpulver]

The out-of-range fix looks correct. Following on from the sum_functions dtype catch above — the same hazard is still live one function down, in Tabulated1D.__call__ itself, and this PR leaves it. That path builds its output with y = np.zeros_like(x) (no dtype=float, unlike the accumulator you just fixed in sum_functions), so it inherits the input array's dtype. Hand it an integer array and y is integer, so every interpolated value — and now the endpoint value your y[idx < 0] / y[idx > len(self.x) - 2] lines assign — is truncated to int on write. Silent: wrong numbers, no error.

All points in-range here, so this is distinct from the out-of-range issue:

f = openmc.data.Tabulated1D([0.0, 10.0], [0.0, 1.0])   # f(x) = x/10
f(5)                       # 0.5                 scalar path, correct
f(np.array([5]))           # array([0])          <- int input, truncated
f(np.array([5.0]))         # array([0.5])        <- float input, correct
f(np.array([2, 5, 8]))     # array([0, 0, 0])    correct: [0.2, 0.5, 0.8]

The output dtype just tracks the input dtype:

input int64    -> [0 0 0]        dtype int64
input float64  -> [0.2 0.5 0.8]  dtype float64

Checked against this branch specifically: the two out-of-range lines don't change it (in-range points never reach them), and for out-of-range integer input the endpoint value they assign is itself truncated.

Same one-word fix you applied to sum_functions:

y = np.zeros_like(x, dtype=float)

In transport this is latent since energies are floats, but the public API accepts integer arrays (the issue's own repro passes them), so it can surface in post-processing/plotting. Happy to push the one-liner plus a small regression test, or it folds cleanly into this PR.

[kawacukennedy]

Thanks @CAOShurong and @dylanpulver for catching the same dtype hazard in Tabulated1D.__call__ — you're exactly right, it was the same latent bug in the caller itself.

Reproduced (the exact case from the review):

f = openmc.data.Tabulated1D([0.0, 10.0], [0.0, 1.0])
f(np.array([5]))    # array([0])   <- int64 truncation, silent wrong numbers
f(np.array([5.0]))  # array([0.5]) <- correct

Fixed (commit 22413cb): changed the output accumulator in Tabulated1D.__call__ from np.zeros_like(x) to np.zeros_like(x, dtype=float), preventing truncation regardless of input dtype.

Added test_tabulated1d_integer_input: reproduces the exact truncation scenario — verifies scalar/array consistency on integer input, output dtype is float64, and values match expected floats.

Full test_function.py suite: 13/13 passed under --noconftest (C++ lib not needed for these pure-Python tests). The sum_functions dtype fix from commit 56a8e5a remains in place — both hazards are now covered.