Tabulated1D array evaluation snaps points within 1e-5 RELATIVE of an endpoint to it (np.isclose default rtol), zeroing cross sections just above thresholds

[shimwell]

Tabulated1D.__call__ snaps any x within 1e-5 relative of the first or last abscissa to that endpoint's y, because np.isclose is called with atol=1e-14 but numpy's default rtol=1e-5 still applies. The scalar path does not do this, so the two disagree, and the array path is the wrong one.

https://github.com/shimwell/endf-python/blob/main/src/endf/function.py#L135-L136

# In some cases, x values might be outside the tabulated region due only
# to precision, so we check if they're close and set them equal if so.
y[np.isclose(x, self.x[0], atol=1e-14)] = self.y[0]
y[np.isclose(x, self.x[-1], atol=1e-14)] = self.y[-1]

The comment says "due only to precision", and atol=1e-14 reads as that intent. But the effective tolerance is 1e-14 + 1e-5 * |x[0]|, which at 6 MeV is 67 eV and at a 5991 eV edge is 0.06 eV. Any genuine grid point inside that window is overwritten with the endpoint value, which at a threshold is zero.

Two confirmed cases

Ba L1 photoionization (photoat-056_Ba_000.endf, MF=23 MT=535), x[0] = 5991.0:

xs(5991.009)                      # 21012.35694510297   scalar, correct
xs(np.array([5991.009]))          # 0.0                 array, snapped
# gap 0.009 eV; window 1e-5 * 5991 = 0.0599 eV

Zero cross section 9 meV above an L1 edge whose value is 21,012 barns.

N15 (n,n3) on its ACE grid, x[0] = 6749151.0, next grid point 6749200.0:

# gap 49 eV; window 1e-5 * 6749151 = 67.5 eV  -> snapped to y[0] = 0

That one point is then missing from the MT 4 inelastic sum, so MT 4 no longer equals the sum of MT 51..91, by 4.61e-6 barns.

Note the MF=3 tabulation for the same reaction is not affected (x[0] = 6747200, gap 2000 eV, window 67 eV), which is why this only appears after NJOY processing. The window scales with the energy, so it bites wherever a grid happens to be finer than 1e-5 relative, and nowhere else. That is why it is rare and easy to miss.

Where it comes from: openmc has it too

These two lines are in openmc as well, at openmc/data/function.py:207-208, identical. endf-python inherited them in 503e4cc, "Initial ENDF functionality ported over from OpenMC", and they are unchanged in paulromano/endf-python main (src/endf/function.py:122-123) and in every fork of it.

Verified on openmc 0.15.3 with the real Ba tabulation, whose first abscissa IS the edge:

import numpy as np, openmc.data
t = openmc.data.Tabulated1D(
    np.array([5991.0, 5991.0, 6167.88862, 6350.0]),
    np.array([0.0, 21012.409, 19989.3069, 19016.02]),
    breakpoints=[4], interpolation=[2])

t(5991.009)                     # 21012.35694510297   scalar
t(np.array([5991.009]))[0]      # 0.0                 array

So this is not a regression against openmc; it is a long-standing behaviour that openmc, endf-python and every downstream consumer share. It matters more here than in openmc because a converter evaluates whole grids at once, which is the array path, whereas openmc's transport reads its own processed tables.

A note on testing it: the window is 1e-5 * |x[0]|, so a synthetic tabulation with a leading point far below the edge will not reproduce it. x[0] has to be at or near the value being probed. My first attempt at this comparison used a leading point at 5000 eV and wrongly concluded openmc was unaffected.

How it was found

Converting ENDF/B-VIII.1 to a different format and diffing against previously published output. Of 33 photoatomic elements compared, Ba was the only one that showed it: only a grid point that lands inside the window is affected, and most grids have none. On the neutron side it showed as N15's MT 4 disagreeing with its own components at one point in 979.

Suggested fix

Pass rtol=0 so the tolerance is the absolute one the comment intends:

y[np.isclose(x, self.x[0], rtol=0, atol=1e-14)] = self.y[0]
y[np.isclose(x, self.x[-1], rtol=0, atol=1e-14)] = self.y[-1]

A regression test worth having: evaluate the same x scalar and array and require agreement, using a tabulation whose second abscissa is within 1e-5 relative of the first.

[shimwell]

this appears to be in the main branch of endf-python

[shimwell]

Correction to the original text: I first wrote that openmc does not have this. That was wrong, and the issue body is now fixed.

My check constructed a tabulation with a leading abscissa at 5000 eV, which makes the window 1e-5 * 5000 = 0.05 eV while the probe sat 991 eV away, so the snapping could never trigger. With the real leading abscissa (x[0] = 5991.0, the edge itself) openmc 0.15.3 returns 0.0 from the array path and 21012.357 from the scalar path, exactly as endf-python does.

So the lines originate in openmc and are shared by openmc, paulromano/endf-python and its forks. Worth raising upstream in openmc as well.

[shimwell]

Widening the check to all 100 ENDF/B-VIII.1 photoatomic elements found two more, and they pin the mechanism quantitatively. In each case the published photoionization cross section is exactly 0.0 at the first grid point above an absorption edge, and in each case that point falls inside the 1e-5 * x[0] window:

element	subshell	edge x[0]	first point above	gap	window	inside
Ba	L1	5991 eV	5991.009290	0.00929 eV	0.05991 eV	yes
Yb	N2	399 eV	399.003783	0.003783 eV	0.00399 eV	yes, by 5%
Ra	M2	4491 eV	4491.006960	0.00696 eV	0.04491 eV	yes

Yb is the sharpest of the three: the gap is 95% of the window, so it is affected and a grid one part in twenty coarser would not have been. That is the signature of a relative tolerance rather than the ulp-scale guard the comment describes, and it is why only 3 of 100 elements show it at all.

97 of the 100 elements are unaffected and match to the last bit.