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.

Addressed

[paulromano]

I took a closer look at the proposed out-of-range behavior, including both the ENDF-6 specification and existing Python callers of Tabulated1D.

ENDF-6 defines the interpolation laws between tabulated points and describes the tabulated interval as the complete region over which the independent variable is defined. It does not prescribe extrapolation below the first point or above the last point. Thus, endpoint clamping and zero extension are both extensions beyond the ENDF-6 format.

The existing Python codebase does rely on the longstanding array behavior of returning zero outside the tabulated interval. In particular, threshold reaction cross sections are evaluated on union grids in ENDF activation and photon-production processing, photoatomic HDF5 export, redundant cross-section construction, and a few other areas. The revised approach is therefore:

1. Preserve the established Python array convention that values outside the tabulated interval are zero.
2. Make scalar evaluation follow the same convention, resolving #4041 without changing existing array-based reaction processing.
3. Document explicitly that zero extension is an OpenMC convention; it is not specified by ENDF-6.
4. Continue returning floating-point arrays for integer-valued inputs, avoiding silent truncation during interpolation.
5. Reject NaN and infinite arguments consistently for scalar and array inputs.
6. Keep sum_functions materializing a Tabulated1D on the union grid, as it did previously. This preserves the HDF5 representation used for derived fission-energy quantities. As before, this materialization cannot represent a discontinuity between a nonzero endpoint and zero extension exactly; addressing that would require a broader cross-language design change.

One related point is that the C++ Tabulated1D class currently clamps to the nearest endpoint. Again, ENDF does not require that behavior, and changing it could affect transport, so I think a cross-language extrapolation policy should be considered separately rather than silently changing the longstanding Python array semantics in this bug fix.