From e396885c4489eff585c2a5dc0f0a02fc2cb93b1d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 13 Oct 2026 11:41:56 -0500 Subject: [PATCH 1/3] Fix MF4 Legendre angular distribution weighting --- CHANGELOG.md | 1 + src/endf/mf4.py | 4 ++ tests/test_mf4.py | 124 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 129 insertions(+) create mode 100644 tests/test_mf4.py diff --git a/CHANGELOG.md b/CHANGELOG.md index b927da3..3a9739d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0. ### Fixed +* Correct Legendre weighting in the MF=4 angular distribution conversion * Fixing MF7/MT2 key name: LHTR → LTHR ([#37](https://github.com/paulromano/endf-python/pull/37)) ## [0.1.12] diff --git a/src/endf/mf4.py b/src/endf/mf4.py index 800dcb3..fce0e70 100644 --- a/src/endf/mf4.py +++ b/src/endf/mf4.py @@ -119,6 +119,8 @@ def from_dict(cls, data: dict) -> AngleDistribution: mu = [] for a_l in data['legendre']['a_l']: coef = np.insert(a_l, 0, 1.0) + # NumPy requires the ENDF angular weighting explicitly. + coef = coef * (2 * np.arange(coef.size) + 1) / 2 mu.append(Legendre(coef)) elif LTT == 2 and LI == 0: energy = data['tabulated']['E'] @@ -129,6 +131,8 @@ def from_dict(cls, data: dict) -> AngleDistribution: mu_leg = [] for a_l in data['legendre']['a_l']: coef = np.insert(a_l, 0, 1.0) + # NumPy requires the ENDF angular weighting explicitly. + coef = coef * (2 * np.arange(coef.size) + 1) / 2 mu_leg.append(Legendre(coef)) # Then get tabulated diff --git a/tests/test_mf4.py b/tests/test_mf4.py new file mode 100644 index 0000000..80b431a --- /dev/null +++ b/tests/test_mf4.py @@ -0,0 +1,124 @@ +# SPDX-FileCopyrightText: 2026 OpenMC contributors +# SPDX-License-Identifier: MIT + +import numpy as np +from numpy.polynomial import Legendre +from numpy.polynomial.legendre import leggauss +import pytest + +from endf.function import Tabulated1D +from endf.mf4 import AngleDistribution + + +def angular_data(ltt, coefficients, energy=None): + if energy is None: + energy = np.arange(1, len(coefficients) + 1, dtype=float) + data = { + 'LTT': ltt, + 'LI': 0, + 'legendre': { + 'E': np.asarray(energy), + 'a_l': [np.array(a_l, dtype=float) for a_l in coefficients], + }, + } + if ltt == 3: + data['tabulated'] = { + 'E': np.array([1.e6, 2.e7]), + 'mu': [ + Tabulated1D([-1., 0., 1.], [0.25, 0.5, 0.75]), + Tabulated1D([-1., 1.], [0.4, 0.6]), + ], + } + return data + + +@pytest.mark.parametrize('ltt', [1, 3]) +@pytest.mark.parametrize( + 'stored, expected_coef, expected_values', + [ + ([], [0.5], [0.5, 0.5, 0.5]), + ([0.2, 0.1], [0.5, 0.3, 0.25], [0.45, 0.375, 1.05]), + ], + ids=['isotropic', 'anisotropic'], +) +def test_legendre_values(ltt, stored, expected_coef, expected_values): + distribution = AngleDistribution.from_dict(angular_data(ltt, [stored])) + density = distribution.mu[0] + + assert isinstance(density, Legendre) + np.testing.assert_allclose(density.coef, expected_coef, rtol=0., atol=1.e-13) + np.testing.assert_allclose( + density(np.array([-1., 0., 1.])), expected_values, rtol=0., atol=1.e-13 + ) + for cosine, expected in zip([-1., 0., 1.], expected_values): + assert density(cosine) == pytest.approx(expected, rel=0., abs=1.e-13) + + +@pytest.mark.parametrize('ltt', [1, 3]) +@pytest.mark.parametrize('stored', [[], [0.2, 0.1]], ids=['isotropic', 'anisotropic']) +def test_legendre_moments(ltt, stored): + distribution = AngleDistribution.from_dict(angular_data(ltt, [stored])) + cosine, weights = leggauss(8) + density = distribution.mu[0](cosine) + + expected_moments = [1., *stored, *([0.] * (2 - len(stored)))] + for order, expected in enumerate(expected_moments): + moment = np.dot(weights, Legendre.basis(order)(cosine) * density) + assert moment == pytest.approx(expected, rel=0., abs=1.e-13) + + +@pytest.mark.parametrize('ltt', [1, 3]) +def test_legendre_multiple_energies_and_input_preservation(ltt): + data = angular_data(ltt, [[], [0.2, 0.1], [-0.1]], [1.e3, 1.e-5, 1.]) + original_energy = data['legendre']['E'].copy() + original_coefficients = [a_l.copy() for a_l in data['legendre']['a_l']] + expected_energy = original_energy + if ltt == 3: + expected_energy = np.hstack((expected_energy, data['tabulated']['E'])) + + first = AngleDistribution.from_dict(data) + second = AngleDistribution.from_dict(data) + np.testing.assert_array_equal(first.energy, expected_energy) + np.testing.assert_array_equal(second.energy, expected_energy) + np.testing.assert_array_equal(data['legendre']['E'], original_energy) + assert len(first.mu) == len(second.mu) == len(expected_energy) + + cosine, weights = leggauss(8) + for i, original in enumerate(original_coefficients): + np.testing.assert_array_equal(data['legendre']['a_l'][i], original) + np.testing.assert_array_equal(first.mu[i].coef, second.mu[i].coef) + expected_moments = [1., *original, *([0.] * (2 - len(original)))] + for order, expected in enumerate(expected_moments): + moment = np.dot( + weights, Legendre.basis(order)(cosine) * first.mu[i](cosine) + ) + assert moment == pytest.approx(expected, rel=0., abs=1.e-13) + + +def test_mixed_tabulated_preservation(): + data = angular_data(3, [[], [0.2, 0.1]]) + tabulated = data['tabulated'] + original_energy = tabulated['E'].copy() + original_values = [(f.x.copy(), f.y.copy()) for f in tabulated['mu']] + + distribution = AngleDistribution.from_dict(data) + + np.testing.assert_array_equal(distribution.energy[2:], original_energy) + np.testing.assert_array_equal(tabulated['E'], original_energy) + for actual, original, (x, y) in zip( + distribution.mu[2:], tabulated['mu'], original_values + ): + assert actual is original + np.testing.assert_array_equal(actual.x, x) + np.testing.assert_array_equal(actual.y, y) + np.testing.assert_array_equal(actual(x), y) + + +def test_pure_tabulated_preservation(): + tabulated = angular_data(3, [[]])['tabulated'] + data = {'LTT': 2, 'LI': 0, 'tabulated': tabulated} + + distribution = AngleDistribution.from_dict(data) + + assert distribution.energy is tabulated['E'] + assert distribution.mu is tabulated['mu'] From 09a8182b799e0f2ec4ff0c47ad9322910f3c7b7c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 13 Oct 2026 12:35:37 -0500 Subject: [PATCH 2/3] Simplify MF4 Legendre regression tests --- src/endf/mf4.py | 6 +- tests/test_mf4.py | 138 +++++++++++++--------------------------------- 2 files changed, 40 insertions(+), 104 deletions(-) diff --git a/src/endf/mf4.py b/src/endf/mf4.py index fce0e70..85a05b8 100644 --- a/src/endf/mf4.py +++ b/src/endf/mf4.py @@ -119,8 +119,7 @@ def from_dict(cls, data: dict) -> AngleDistribution: mu = [] for a_l in data['legendre']['a_l']: coef = np.insert(a_l, 0, 1.0) - # NumPy requires the ENDF angular weighting explicitly. - coef = coef * (2 * np.arange(coef.size) + 1) / 2 + coef *= (2 * np.arange(coef.size) + 1) / 2 mu.append(Legendre(coef)) elif LTT == 2 and LI == 0: energy = data['tabulated']['E'] @@ -131,8 +130,7 @@ def from_dict(cls, data: dict) -> AngleDistribution: mu_leg = [] for a_l in data['legendre']['a_l']: coef = np.insert(a_l, 0, 1.0) - # NumPy requires the ENDF angular weighting explicitly. - coef = coef * (2 * np.arange(coef.size) + 1) / 2 + coef *= (2 * np.arange(coef.size) + 1) / 2 mu_leg.append(Legendre(coef)) # Then get tabulated diff --git a/tests/test_mf4.py b/tests/test_mf4.py index 80b431a..542d8cb 100644 --- a/tests/test_mf4.py +++ b/tests/test_mf4.py @@ -2,123 +2,61 @@ # SPDX-License-Identifier: MIT import numpy as np -from numpy.polynomial import Legendre from numpy.polynomial.legendre import leggauss import pytest - from endf.function import Tabulated1D from endf.mf4 import AngleDistribution -def angular_data(ltt, coefficients, energy=None): - if energy is None: - energy = np.arange(1, len(coefficients) + 1, dtype=float) +@pytest.mark.parametrize('ltt', [1, 3]) +def test_legendre_conversion(ltt): + coefficients = [np.array([]), np.array([0.2, 0.1]), np.array([-0.1])] + energy = np.array([1.e3, 1.e-5, 1.]) + tabulated = Tabulated1D([-1., 0., 1.], [0.25, 0.5, 0.75]) data = { 'LTT': ltt, 'LI': 0, - 'legendre': { - 'E': np.asarray(energy), - 'a_l': [np.array(a_l, dtype=float) for a_l in coefficients], - }, + 'legendre': {'E': energy, 'a_l': coefficients}, + 'tabulated': {'E': np.array([1.e6]), 'mu': [tabulated]}, } - if ltt == 3: - data['tabulated'] = { - 'E': np.array([1.e6, 2.e7]), - 'mu': [ - Tabulated1D([-1., 0., 1.], [0.25, 0.5, 0.75]), - Tabulated1D([-1., 1.], [0.4, 0.6]), - ], - } - return data - - -@pytest.mark.parametrize('ltt', [1, 3]) -@pytest.mark.parametrize( - 'stored, expected_coef, expected_values', - [ - ([], [0.5], [0.5, 0.5, 0.5]), - ([0.2, 0.1], [0.5, 0.3, 0.25], [0.45, 0.375, 1.05]), - ], - ids=['isotropic', 'anisotropic'], -) -def test_legendre_values(ltt, stored, expected_coef, expected_values): - distribution = AngleDistribution.from_dict(angular_data(ltt, [stored])) - density = distribution.mu[0] - - assert isinstance(density, Legendre) - np.testing.assert_allclose(density.coef, expected_coef, rtol=0., atol=1.e-13) - np.testing.assert_allclose( - density(np.array([-1., 0., 1.])), expected_values, rtol=0., atol=1.e-13 - ) - for cosine, expected in zip([-1., 0., 1.], expected_values): - assert density(cosine) == pytest.approx(expected, rel=0., abs=1.e-13) + original_coefficients = [a_l.copy() for a_l in coefficients] + distribution = AngleDistribution.from_dict(data) + repeated = AngleDistribution.from_dict(data) + expected_energy = [1.e3, 1.e-5, 1., 1.e6] if ltt == 3 else [1.e3, 1.e-5, 1.] + np.testing.assert_array_equal(distribution.energy, expected_energy) + np.testing.assert_array_equal(repeated.energy, expected_energy) + np.testing.assert_array_equal(energy, [1.e3, 1.e-5, 1.]) + assert len(distribution.mu) == len(repeated.mu) == len(expected_energy) -@pytest.mark.parametrize('ltt', [1, 3]) -@pytest.mark.parametrize('stored', [[], [0.2, 0.1]], ids=['isotropic', 'anisotropic']) -def test_legendre_moments(ltt, stored): - distribution = AngleDistribution.from_dict(angular_data(ltt, [stored])) + expected_coef = [[0.5], [0.5, 0.3, 0.25], [0.5, -0.15]] + expected_values = [[0.5, 0.5, 0.5], [0.45, 0.375, 1.05], [0.65, 0.5, 0.35]] + expected_moments = [[1., 0., 0.], [1., 0.2, 0.1], [1., -0.1, 0.]] cosine, weights = leggauss(8) - density = distribution.mu[0](cosine) - - expected_moments = [1., *stored, *([0.] * (2 - len(stored)))] - for order, expected in enumerate(expected_moments): - moment = np.dot(weights, Legendre.basis(order)(cosine) * density) - assert moment == pytest.approx(expected, rel=0., abs=1.e-13) + basis = np.array([np.ones(8), cosine, (3 * cosine**2 - 1) / 2]) + for i, density in enumerate(distribution.mu[:3]): + assert density.coef == pytest.approx(expected_coef[i], rel=0., abs=1.e-13) + assert density([-1., 0., 1.]) == pytest.approx(expected_values[i], rel=0., abs=1.e-13) + assert density(0.) == pytest.approx(expected_values[i][1], rel=0., abs=1.e-13) + moments = basis @ (weights * density(cosine)) + assert moments == pytest.approx(expected_moments[i], rel=0., abs=1.e-13) + np.testing.assert_array_equal(coefficients[i], original_coefficients[i]) + np.testing.assert_array_equal(density.coef, repeated.mu[i].coef) - -@pytest.mark.parametrize('ltt', [1, 3]) -def test_legendre_multiple_energies_and_input_preservation(ltt): - data = angular_data(ltt, [[], [0.2, 0.1], [-0.1]], [1.e3, 1.e-5, 1.]) - original_energy = data['legendre']['E'].copy() - original_coefficients = [a_l.copy() for a_l in data['legendre']['a_l']] - expected_energy = original_energy if ltt == 3: - expected_energy = np.hstack((expected_energy, data['tabulated']['E'])) - - first = AngleDistribution.from_dict(data) - second = AngleDistribution.from_dict(data) - np.testing.assert_array_equal(first.energy, expected_energy) - np.testing.assert_array_equal(second.energy, expected_energy) - np.testing.assert_array_equal(data['legendre']['E'], original_energy) - assert len(first.mu) == len(second.mu) == len(expected_energy) - - cosine, weights = leggauss(8) - for i, original in enumerate(original_coefficients): - np.testing.assert_array_equal(data['legendre']['a_l'][i], original) - np.testing.assert_array_equal(first.mu[i].coef, second.mu[i].coef) - expected_moments = [1., *original, *([0.] * (2 - len(original)))] - for order, expected in enumerate(expected_moments): - moment = np.dot( - weights, Legendre.basis(order)(cosine) * first.mu[i](cosine) - ) - assert moment == pytest.approx(expected, rel=0., abs=1.e-13) - - -def test_mixed_tabulated_preservation(): - data = angular_data(3, [[], [0.2, 0.1]]) - tabulated = data['tabulated'] - original_energy = tabulated['E'].copy() - original_values = [(f.x.copy(), f.y.copy()) for f in tabulated['mu']] - - distribution = AngleDistribution.from_dict(data) - - np.testing.assert_array_equal(distribution.energy[2:], original_energy) - np.testing.assert_array_equal(tabulated['E'], original_energy) - for actual, original, (x, y) in zip( - distribution.mu[2:], tabulated['mu'], original_values - ): - assert actual is original - np.testing.assert_array_equal(actual.x, x) - np.testing.assert_array_equal(actual.y, y) - np.testing.assert_array_equal(actual(x), y) + assert distribution.mu[-1] is repeated.mu[-1] is tabulated + np.testing.assert_array_equal(data['tabulated']['E'], [1.e6]) + np.testing.assert_array_equal(tabulated.x, [-1., 0., 1.]) + np.testing.assert_array_equal(tabulated.y, [0.25, 0.5, 0.75]) + assert tabulated([-1., 0., 1.]) == pytest.approx([0.25, 0.5, 0.75]) -def test_pure_tabulated_preservation(): - tabulated = angular_data(3, [[]])['tabulated'] - data = {'LTT': 2, 'LI': 0, 'tabulated': tabulated} +def test_tabulated_conversion(): + energy = np.array([1.e6]) + mu = [Tabulated1D([-1., 1.], [0.4, 0.6])] + data = {'LTT': 2, 'LI': 0, 'tabulated': {'E': energy, 'mu': mu}} distribution = AngleDistribution.from_dict(data) - assert distribution.energy is tabulated['E'] - assert distribution.mu is tabulated['mu'] + assert distribution.energy is energy + assert distribution.mu is mu From be8b6b60d950e3951eaa70f45ec18d4aef81f9ea Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 13 Oct 2026 14:10:54 -0500 Subject: [PATCH 3/3] Update CHANGELOG --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3a9739d..8f22c96 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,7 +13,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0. ### Fixed -* Correct Legendre weighting in the MF=4 angular distribution conversion +* Correct Legendre weighting in the MF=4 angular distribution conversion ([#46](https://github.com/paulromano/endf-python/pull/46)) * Fixing MF7/MT2 key name: LHTR → LTHR ([#37](https://github.com/paulromano/endf-python/pull/37)) ## [0.1.12]