from math import pi
from pathlib import Path

import numpy as np
import pytest
import openmc
import openmc.stats
from openmc.stats.univariate import _INTERPOLATION_SCHEMES, DecaySpectrum
from scipy.integrate import trapezoid

from tests.unit_tests import assert_sample_mean


@pytest.mark.flaky(reruns=1)
def test_discrete():
    x = [0.0, 1.0, 10.0]
    p = [0.3, 0.2, 0.5]
    d = openmc.stats.Discrete(x, p)
    elem = d.to_xml_element('distribution')

    d = openmc.stats.Discrete.from_xml_element(elem)
    np.testing.assert_array_equal(d.x, x)
    np.testing.assert_array_equal(d.p, p)
    assert len(d) == len(x)

    d = openmc.stats.Univariate.from_xml_element(elem)
    assert isinstance(d, openmc.stats.Discrete)

    # Single point
    d2 = openmc.stats.Discrete(1e6, 1.0)
    assert d2.x == [1e6]
    assert d2.p == [1.0]
    assert len(d2) == 1

    vals = np.array([1.0, 2.0, 3.0])
    probs = np.array([0.1, 0.7, 0.2])

    exp_mean = (vals * probs).sum()

    d3 = openmc.stats.Discrete(vals, probs)

    # sample discrete distribution and check that the mean of the samples is
    # within 4 std. dev. of the expected mean
    n_samples = 1_000_000
    samples, weights = d3.sample(n_samples)
    assert_sample_mean(samples, exp_mean)
    assert np.all(weights == 1.0)

    # Test biased distribution
    d3.bias = np.array([0.2, 0.1, 0.7])
    bias_elem = d3.to_xml_element('distribution')
    d4 = openmc.stats.Univariate.from_xml_element(bias_elem)
    np.testing.assert_array_equal(d4.bias, [0.2, 0.1, 0.7])
    samples, weights = d4.sample(n_samples)
    weighted_sample = samples * weights
    assert_sample_mean(weighted_sample, exp_mean)
    assert np.all(weights != 1.0)


def test_delta_function():
    d = openmc.stats.delta_function(14.1e6)
    assert isinstance(d, openmc.stats.Discrete)
    np.testing.assert_array_equal(d.x, [14.1e6])
    np.testing.assert_array_equal(d.p, [1.0])


def test_merge_discrete():
    x1 = [0.0, 1.0, 10.0]
    p1 = [0.3, 0.2, 0.5]
    d1 = openmc.stats.Discrete(x1, p1)

    x2 = [0.5, 1.0, 5.0]
    p2 = [0.4, 0.5, 0.1]
    d2 = openmc.stats.Discrete(x2, p2)

    # Merged distribution should have x values sorted and probabilities
    # appropriately combined. Duplicate x values should appear once.
    merged = openmc.stats.Discrete.merge([d1, d2], [0.6, 0.4])
    assert merged.x == pytest.approx([0.0, 0.5, 1.0, 5.0, 10.0])
    assert merged.p == pytest.approx(
        [0.6*0.3, 0.4*0.4, 0.6*0.2 + 0.4*0.5, 0.4*0.1, 0.6*0.5])
    assert merged.integral() == pytest.approx(1.0)

    # Probabilities add up but are not normalized
    d1 = openmc.stats.Discrete([3.0], [1.0])
    triple = openmc.stats.Discrete.merge([d1, d1, d1], [1.0, 2.0, 3.0])
    assert triple.x == pytest.approx([3.0])
    assert triple.p == pytest.approx([6.0])
    assert triple.integral() == pytest.approx(6.0)


def test_merge_discrete_with_bias():
    # Two discrete distributions with different biases
    d1 = openmc.stats.Discrete([1.0, 2.0], [0.5, 0.5])
    d2 = openmc.stats.Discrete([2.0, 3.0], [0.3, 0.7], bias=[0.1, 0.9])

    merged = openmc.stats.Discrete.merge([d1, d2], [0.6, 0.4])
    exp_mean = 0.6 * d1.mean() + 0.4 * d2.mean()

    # Verify merged distribution has correct x values
    assert set(merged.x) == {1.0, 2.0, 3.0}

    # Bias should not be changed in original distributions
    assert d1.bias is None
    assert np.all(d2.bias == [0.1, 0.9])

    # Sample and verify bias is applied correctly
    samples, weights = merged.sample(10_000)

    # Verify weighted mean matches expected unbiased mean
    assert_sample_mean(samples*weights, exp_mean)


def test_clip_discrete():
    # Create discrete distribution with two points that are not important, one
    # because the x value is very small, and one because the p value is very
    # small
    d = openmc.stats.Discrete([1e-8, 1.0, 2.0, 1000.0], [3.0, 2.0, 5.0, 1e-12])

    # Clipping the distribution should result in two points
    d_clip = d.clip(1e-6)
    assert d_clip.x.size == 2
    assert d_clip.p.size == 2

    # Make sure inplace returns same object
    d_same = d.clip(1e-6, inplace=True)
    assert d_same is d

    with pytest.raises(ValueError):
        d.clip(-1.)

    with pytest.raises(ValueError):
        d.clip(5)


@pytest.mark.flaky(reruns=1)
def test_uniform():
    a, b = 10.0, 20.0
    d = openmc.stats.Uniform(a, b)
    elem = d.to_xml_element('distribution')

    d = openmc.stats.Uniform.from_xml_element(elem)
    assert d.a == a
    assert d.b == b
    assert len(d) == 2

    t = d.to_tabular()
    np.testing.assert_array_equal(t.x, [a, b])
    np.testing.assert_array_equal(t.p, [1/(b-a), 1/(b-a)])
    assert t.interpolation == 'histogram'

    # Sample distribution and check that the mean of the samples is within 4
    # std. dev. of the expected mean
    exp_mean = 0.5 * (a + b)
    n_samples = 1_000_000
    samples, weights = d.sample(n_samples)
    assert_sample_mean(samples, exp_mean)
    assert np.all(weights == 1.0)

    # Test biased distribution
    d.bias = openmc.stats.PowerLaw(a, b, 2)
    bias_elem = d.to_xml_element('distribution')
    d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
    assert isinstance (d2.bias, openmc.stats.PowerLaw)
    samples, weights = d2.sample(n_samples)
    weighted_sample = samples * weights
    assert_sample_mean(weighted_sample, exp_mean)
    assert np.all(weights != 1.0)


@pytest.mark.flaky(reruns=1)
def test_powerlaw():
    a, b, n = 10.0, 100.0, 2.0
    d = openmc.stats.PowerLaw(a, b, n)
    elem = d.to_xml_element('distribution')

    d = openmc.stats.PowerLaw.from_xml_element(elem)
    assert d.a == a
    assert d.b == b
    assert d.n == n
    assert len(d) == 3

    # Determine mean of distribution
    exp_mean = (n+1)*(b**(n+2) - a**(n+2))/((n+2)*(b**(n+1) - a**(n+1)))

    # sample power law distribution and check that the mean of the samples is
    # within 4 std. dev. of the expected mean
    n_samples = 1_000_000
    samples, weights = d.sample(n_samples)
    assert_sample_mean(samples, exp_mean)
    assert np.all(weights == 1.0)

    # Test biased distribution
    d.bias = openmc.stats.Uniform(a, b)
    bias_elem = d.to_xml_element('distribution')
    d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
    assert isinstance (d2.bias, openmc.stats.Uniform)
    samples, weights = d2.sample(n_samples)
    weighted_sample = samples * weights
    assert_sample_mean(weighted_sample, exp_mean)
    assert np.all(weights != 1.0)


@pytest.mark.flaky(reruns=1)
def test_maxwell():
    theta = 1.2895e6
    d = openmc.stats.Maxwell(theta)
    elem = d.to_xml_element('distribution')

    d = openmc.stats.Maxwell.from_xml_element(elem)
    assert d.theta == theta
    assert len(d) == 1

    exp_mean = 3/2 * theta

    # sample maxwell distribution and check that the mean of the samples is
    # within 4 std. dev. of the expected mean
    n_samples = 1_000_000
    samples, weights = d.sample(n_samples)
    assert_sample_mean(samples, exp_mean)
    assert np.all(weights == 1.0)

    # A second sample starting from a different seed
    samples_2, weights_2 = d.sample(n_samples)
    assert_sample_mean(samples_2, exp_mean)
    assert samples_2.mean() != samples.mean()
    assert np.all(weights_2 == 1)

    # Test biased distribution
    d.bias = openmc.stats.Maxwell((theta * 1.1))
    bias_elem = d.to_xml_element('distribution')
    d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
    assert isinstance (d2.bias, openmc.stats.Maxwell)
    samples, weights = d2.sample(n_samples)
    weighted_sample = samples * weights
    assert_sample_mean(weighted_sample, exp_mean)
    assert np.all(weights != 1.0)


@pytest.mark.flaky(reruns=1)
def test_watt():
    a, b = 0.965e6, 2.29e-6
    d = openmc.stats.Watt(a, b)
    elem = d.to_xml_element('distribution')

    d = openmc.stats.Watt.from_xml_element(elem)
    assert d.a == a
    assert d.b == b
    assert len(d) == 2

    # mean value form adapted from
    # "Prompt-fission-neutron average energy for 238U(n, f ) from
    # threshold to 200 MeV" Ethvignot et. al.
    # https://doi.org/10.1016/j.physletb.2003.09.048
    exp_mean = 3/2 * a + a**2 * b / 4

    # sample Watt distribution and check that the mean of the samples is within
    # 4 std. dev. of the expected mean
    n_samples = 1_000_000
    samples, weights = d.sample(n_samples)
    assert_sample_mean(samples, exp_mean)
    assert np.all(weights == 1.0)

    # Test biased distribution with 5 percent higher T_e
    d.bias = openmc.stats.Watt(a*1.05, b)
    bias_elem = d.to_xml_element('distribution')
    d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
    assert isinstance (d2.bias, openmc.stats.Watt)
    samples, weights = d2.sample(n_samples)
    weighted_sample = samples * weights
    assert_sample_mean(weighted_sample, exp_mean)
    assert np.all(weights != 1.0)


@pytest.mark.flaky(reruns=1)
def test_tabular():
    # test linear-linear sampling
    x = np.array([0.001, 5.0, 7.0, 10.0])
    p = np.array([10.0, 20.0, 5.0, 6.0])
    d = openmc.stats.Tabular(x, p, 'linear-linear')
    n_samples = 100_000
    samples, weights = d.sample(n_samples)
    assert_sample_mean(samples, d.mean())
    assert np.all(weights == 1.0)

    for scheme in _INTERPOLATION_SCHEMES:
        # test sampling
        d = openmc.stats.Tabular(x, p, scheme)
        n_samples = 100_000
        samples = d.sample(n_samples)[0]
        assert_sample_mean(samples, d.mean())

    # test histogram sampling
    d = openmc.stats.Tabular(x, p, interpolation='histogram')
    samples, weights = d.sample(n_samples)
    assert_sample_mean(samples, d.mean())
    assert np.all(weights == 1.0)

    # Multiplying the probabilities should preserve the mean but change the integral
    d2 = openmc.stats.Tabular(x, p*2, interpolation='histogram')
    assert d2.mean() == pytest.approx(d.mean())
    assert d2.integral() == pytest.approx(2.0*d.integral())

    # Normalizing should result in an integral of 1
    d.normalize()
    assert d.integral() == pytest.approx(1.0)

    # ensure that passing a set of probabilities shorter than x works
    # for histogram interpolation
    d = openmc.stats.Tabular(x, p[:-1], interpolation='histogram')
    d.cdf()
    d.mean()
    assert_sample_mean(d.sample(n_samples)[0], d.mean())

    # passing a shorter probability set should raise an error for linear-linear
    with pytest.raises(ValueError):
        d = openmc.stats.Tabular(x, p[:-1], interpolation='linear-linear')
        d.cdf()

    # Use probabilities of correct length for linear-linear interpolation and
    # call the CDF method
    d = openmc.stats.Tabular(x, p, interpolation='linear-linear')
    d.cdf()

    # Test biased distribution
    d.bias = openmc.stats.Uniform(x[0], x[-1])
    bias_elem = d.to_xml_element('distribution')
    d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
    assert isinstance (d2.bias, openmc.stats.Uniform)
    samples, weights = d2.sample(n_samples)
    weighted_sample = samples * weights
    assert_sample_mean(weighted_sample, d2.mean())
    assert np.all(weights != 1.0)


def test_tabular_from_xml():
    x = np.array([0.0, 5.0, 7.0, 10.0])
    p = np.array([10.0, 20.0, 5.0, 6.0])
    d = openmc.stats.Tabular(x, p, 'linear-linear')
    elem = d.to_xml_element('distribution')

    d = openmc.stats.Tabular.from_xml_element(elem)
    assert all(d.x == x)
    assert all(d.p == p)
    assert d.interpolation == 'linear-linear'
    assert len(d) == len(x)

    # Make sure XML roundtrip works with len(x) == len(p) + 1
    x = np.array([0.0, 5.0, 7.0, 10.0])
    p = np.array([10.0, 20.0, 5.0])
    d = openmc.stats.Tabular(x, p, 'histogram')
    elem = d.to_xml_element('distribution')
    d = openmc.stats.Tabular.from_xml_element(elem)
    assert all(d.x == x)
    assert all(d.p == p)


def test_legendre():
    # Pu239 elastic scattering at 100 keV
    coeffs = [1.000e+0, 1.536e-1, 1.772e-2, 5.945e-4, 3.497e-5, 1.881e-5]
    d = openmc.stats.Legendre(coeffs)
    assert d.coefficients == pytest.approx(coeffs)
    assert len(d) == len(coeffs)

    # Integrating distribution should yield one
    mu = np.linspace(-1., 1., 1000)
    assert trapezoid(d(mu), mu) == pytest.approx(1.0, rel=1e-4)

    with pytest.raises(NotImplementedError):
        d.to_xml_element('distribution')


@pytest.mark.flaky(reruns=1)
def test_mixture():
    d1 = openmc.stats.Uniform(0, 5)
    d2 = openmc.stats.Uniform(3, 7)
    p = [0.5, 0.5]
    mix = openmc.stats.Mixture(p, [d1, d2])
    np.testing.assert_allclose(mix.probability, p)
    assert mix.distribution == [d1, d2]
    assert len(mix) == 4

    # Sample and make sure sample mean is close to expected mean
    n_samples = 1_000_000
    samples, weights = mix.sample(n_samples)
    assert_sample_mean(samples, (2.5 + 5.0)/2)
    assert np.all(weights == 1.0)

    elem = mix.to_xml_element('distribution')

    d = openmc.stats.Mixture.from_xml_element(elem)
    np.testing.assert_allclose(d.probability, p)
    assert d.distribution == [d1, d2]
    assert len(d) == 4

    # Test biased sub-distribution
    d.distribution[0].bias = openmc.stats.PowerLaw(0, 5, 2)
    bias_elem = d.to_xml_element('distribution')
    d3 = openmc.stats.Univariate.from_xml_element(bias_elem)
    assert isinstance (d3.distribution[0].bias, openmc.stats.PowerLaw)
    samples, weights = d3.sample(n_samples)
    weighted_sample = samples * weights
    assert_sample_mean(weighted_sample, (2.5 + 5.0)/2)

    # Test biased meta-probability
    d.distribution[0].bias = None
    d.bias = [0.25, 0.75]
    bias_elem_2 = d.to_xml_element('distribution')
    d4 = openmc.stats.Univariate.from_xml_element(bias_elem_2)
    assert isinstance (d4.bias, np.ndarray)
    samples, weights = d4.sample(n_samples)
    weighted_sample = samples * weights
    assert_sample_mean(weighted_sample, (2.5 + 5.0)/2)
    assert np.all(weights != 1.0)


def test_mixture_clip():
    # Create mixture distribution containing a discrete distribution with two
    # points that are not important, one because the x value is very small, and
    # one because the p value is very small
    d1 = openmc.stats.Discrete([1e-8, 1.0, 2.0, 1000.0], [3.0, 2.0, 5.0, 1e-12])
    d2 = openmc.stats.Uniform(0, 5)
    mix = openmc.stats.Mixture([0.5, 0.5], [d1, d2])

    # Clipping should reduce the contained discrete distribution to 2 points
    mix_clip = mix.clip(1e-6)
    assert mix_clip.distribution[0].x.size == 2
    assert mix_clip.distribution[0].p.size == 2

    # Make sure inplace returns same object
    mix_same = mix.clip(1e-6, inplace=True)
    assert mix_same is mix

    # Make sure clip removes low probability distributions
    d_small = openmc.stats.Uniform(0., 1.)
    d_large = openmc.stats.Uniform(2., 5.)
    mix = openmc.stats.Mixture([1e-10, 1.0], [d_small, d_large])
    mix_clip = mix.clip(1e-3)
    assert mix_clip.distribution == [d_large]

    # Make sure warning is raised if tolerance is exceeded
    d1 = openmc.stats.Discrete([1.0, 1.001], [1.0, 0.7e-6])
    d2 = openmc.stats.Tabular([0.0, 1.0], [0.7e-6], interpolation='histogram')
    mix = openmc.stats.Mixture([1.0, 1.0], [d1, d2])
    with pytest.warns(UserWarning):
        mix_clip = mix.clip(1e-6)

    # Make sure warning is raised if a biased Discrete is clipped
    d3 = openmc.stats.Discrete([1.0, 1.001], [1.0, 0.7e-8])
    d3.bias = [0.9, 0.1]
    mix = openmc.stats.Mixture([1.0, 1.0], [d3, d2])
    with pytest.raises(RuntimeError):
        mix_clip = mix.clip(1e-6)


def test_polar_azimuthal():
    # default polar-azimuthal should be uniform in mu and phi
    d = openmc.stats.PolarAzimuthal()
    assert isinstance(d.mu, openmc.stats.Uniform)
    assert d.mu.a == -1.
    assert d.mu.b == 1.
    assert isinstance(d.phi, openmc.stats.Uniform)
    assert d.phi.a == 0.
    assert d.phi.b == 2*pi

    mu = openmc.stats.Discrete(1., 1.)
    phi = openmc.stats.Discrete(0., 1.)
    d = openmc.stats.PolarAzimuthal(mu, phi)
    assert d.mu == mu
    assert d.phi == phi

    elem = d.to_xml_element()
    assert elem.tag == 'angle'
    assert elem.attrib['type'] == 'mu-phi'
    assert elem.find('mu') is not None
    assert elem.find('phi') is not None

    d = openmc.stats.PolarAzimuthal.from_xml_element(elem)
    assert d.mu.x == [1.]
    assert d.mu.p == [1.]
    assert d.phi.x == [0.]
    assert d.phi.p == [1.]

    d = openmc.stats.UnitSphere.from_xml_element(elem)
    assert isinstance(d, openmc.stats.PolarAzimuthal)


def test_isotropic():
    d = openmc.stats.Isotropic()
    mu = openmc.stats.Uniform(-1.0, 1.0)
    phi = openmc.stats.PowerLaw(0., 2*np.pi, 2)
    d2 = openmc.stats.PolarAzimuthal(mu, phi)
    d.bias = d2
    elem = d.to_xml_element()
    assert elem.tag == 'angle'
    assert elem.attrib['type'] == 'isotropic'

    d = openmc.stats.Isotropic.from_xml_element(elem)
    assert isinstance(d, openmc.stats.Isotropic)
    assert isinstance(d.bias, openmc.stats.PolarAzimuthal)


def test_monodirectional():
    d = openmc.stats.Monodirectional((1., 0., 0.))
    elem = d.to_xml_element()
    assert elem.tag == 'angle'
    assert elem.attrib['type'] == 'monodirectional'

    d = openmc.stats.Monodirectional.from_xml_element(elem)
    assert d.reference_uvw == pytest.approx((1., 0., 0.))


def test_cartesian():
    x = openmc.stats.Uniform(-10., 10.)
    y = openmc.stats.Uniform(-10., 10.)
    z = openmc.stats.Uniform(0., 20.)
    z.bias = openmc.stats.PowerLaw(0., 20., 3)
    d = openmc.stats.CartesianIndependent(x, y, z)

    elem = d.to_xml_element()
    assert elem.tag == 'space'
    assert elem.attrib['type'] == 'cartesian'
    assert elem.find('x') is not None
    assert elem.find('y') is not None

    d = openmc.stats.CartesianIndependent.from_xml_element(elem)
    assert d.x == x
    assert d.y == y
    assert d.z == z

    d = openmc.stats.Spatial.from_xml_element(elem)
    assert isinstance(d, openmc.stats.CartesianIndependent)
    assert isinstance (d.z.bias, openmc.stats.PowerLaw)


def test_box():
    lower_left = (-10., -10., -10.)
    upper_right = (10., 10., 10.)
    d = openmc.stats.Box(lower_left, upper_right)

    elem = d.to_xml_element()
    assert elem.tag == 'space'
    assert elem.attrib['type'] == 'box'
    assert elem.find('parameters') is not None

    d = openmc.stats.Box.from_xml_element(elem)
    assert d.lower_left == pytest.approx(lower_left)
    assert d.upper_right == pytest.approx(upper_right)


def test_point():
    p = (-4., 2., 10.)
    d = openmc.stats.Point(p)

    elem = d.to_xml_element()
    assert elem.tag == 'space'
    assert elem.attrib['type'] == 'point'
    assert elem.find('parameters') is not None

    d = openmc.stats.Point.from_xml_element(elem)
    assert d.xyz == pytest.approx(p)


def test_spherical_uniform():
    r_outer = 2.0
    r_inner = 1.0
    thetas = (0.0, pi/2)
    phis = (0.0, pi)
    origin = (0.0, 1.0, 2.0)

    sph_indep_function = openmc.stats.spherical_uniform(r_outer,
                                                        r_inner,
                                                        thetas,
                                                        phis,
                                                        origin)

    assert isinstance(sph_indep_function, openmc.stats.SphericalIndependent)


def test_cylindrical_uniform():
    r_outer = 2.0
    r_inner = 1.0
    height = 1.0
    phis = (0.0, pi)
    origin = (0.0, 1.0, 2.0)

    dist = openmc.stats.cylindrical_uniform(r_outer, height, r_inner, phis,
                                            origin=origin)

    assert isinstance(dist, openmc.stats.CylindricalIndependent)

    # Check r distribution (PowerLaw with exponent 1 for uniform area sampling)
    assert isinstance(dist.r, openmc.stats.PowerLaw)
    assert dist.r.a == pytest.approx(r_inner)
    assert dist.r.b == pytest.approx(r_outer)
    assert dist.r.n == pytest.approx(1.0)

    # Check phi distribution
    assert isinstance(dist.phi, openmc.stats.Uniform)
    assert dist.phi.a == pytest.approx(phis[0])
    assert dist.phi.b == pytest.approx(phis[1])

    # Check z distribution (centered on origin along z_dir)
    assert isinstance(dist.z, openmc.stats.Uniform)
    assert dist.z.a == pytest.approx(-height / 2)
    assert dist.z.b == pytest.approx(height / 2)

    # Check origin and default directions
    np.testing.assert_allclose(dist.origin, origin)
    np.testing.assert_allclose(dist.r_dir, [1., 0., 0.])
    np.testing.assert_allclose(dist.z_dir, [0., 0., 1.])

    # XML round-trip preserves all parameters
    elem = dist.to_xml_element()
    dist2 = openmc.stats.CylindricalIndependent.from_xml_element(elem)
    np.testing.assert_allclose(dist2.origin, origin)
    np.testing.assert_allclose(dist2.r_dir, dist.r_dir)
    np.testing.assert_allclose(dist2.z_dir, dist.z_dir)


def test_cylindrical_uniform_tilted():
    # Test with non-default axis orientation (y-axis as cylinder axis)
    dist = openmc.stats.cylindrical_uniform(
        r_outer=3.0, height=2.0, r_dir=(1., 0., 0.), z_dir=(0., 1., 0.)
    )
    np.testing.assert_allclose(dist.z_dir, [0., 1., 0.])
    np.testing.assert_allclose(dist.r_dir, [1., 0., 0.])

    # XML round-trip preserves tilted directions
    elem = dist.to_xml_element()
    dist2 = openmc.stats.CylindricalIndependent.from_xml_element(elem)
    np.testing.assert_allclose(dist2.z_dir, dist.z_dir)
    np.testing.assert_allclose(dist2.r_dir, dist.r_dir)


def test_cylindrical_uniform_ring():
    # height=0 should produce a flat ring (delta function at z=0)
    r_outer = 2.0
    r_inner = 1.0
    phis = (0.0, pi)
    origin = (0.0, 1.0, 2.0)

    dist = openmc.stats.cylindrical_uniform(r_outer, 0.0, r_inner, phis,
                                            origin=origin)

    assert isinstance(dist, openmc.stats.CylindricalIndependent)

    # Check r distribution
    assert isinstance(dist.r, openmc.stats.PowerLaw)
    assert dist.r.a == pytest.approx(r_inner)
    assert dist.r.b == pytest.approx(r_outer)
    assert dist.r.n == pytest.approx(1.0)

    # Check phi distribution
    assert isinstance(dist.phi, openmc.stats.Uniform)
    assert dist.phi.a == pytest.approx(phis[0])
    assert dist.phi.b == pytest.approx(phis[1])

    # z distribution must be a delta function at 0.0 (local frame)
    assert isinstance(dist.z, openmc.stats.Discrete)
    assert dist.z.x[0] == pytest.approx(0.0)

    # XML round-trip
    elem = dist.to_xml_element()
    dist2 = openmc.stats.CylindricalIndependent.from_xml_element(elem)
    np.testing.assert_allclose(dist2.origin, origin)
    np.testing.assert_allclose(dist2.r_dir, dist.r_dir)
    np.testing.assert_allclose(dist2.z_dir, dist.z_dir)


@pytest.mark.flaky(reruns=1)
def test_normal():
    mean = 10.0
    std_dev = 2.0
    d = openmc.stats.Normal(mean,std_dev)

    elem = d.to_xml_element('distribution')
    assert elem.attrib['type'] == 'normal'

    d = openmc.stats.Normal.from_xml_element(elem)
    assert d.mean_value == pytest.approx(mean)
    assert d.std_dev == pytest.approx(std_dev)
    assert len(d) == 2

    # sample normal distribution
    n_samples = 100_000
    samples, weights = d.sample(n_samples)
    assert_sample_mean(samples, mean)
    assert np.all(weights == 1.0)

    # Test biased distribution
    d.bias = openmc.stats.Normal(10.0, 4.0)
    bias_elem = d.to_xml_element('distribution')
    d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
    assert isinstance (d2.bias, openmc.stats.Normal)
    samples, weights = d2.sample(n_samples)
    weighted_sample = samples * weights
    assert_sample_mean(weighted_sample, mean)
    assert np.all(weights != 1.0)


@pytest.mark.flaky(reruns=1)
def test_normal_truncated():
    mean = 10.0
    std_dev = 2.0
    lower = 6.0
    upper = 14.0

    d = openmc.stats.Normal(mean, std_dev, lower, upper)

    # Check attributes
    assert d.mean_value == pytest.approx(mean)
    assert d.std_dev == pytest.approx(std_dev)
    assert d.lower == pytest.approx(lower)
    assert d.upper == pytest.approx(upper)
    assert len(d) == 4
    assert d.support == (lower, upper)

    # Test XML round-trip
    elem = d.to_xml_element('distribution')
    assert elem.attrib['type'] == 'normal'
    params = elem.attrib['parameters'].split()
    assert len(params) == 4

    d2 = openmc.stats.Normal.from_xml_element(elem)
    assert d2.mean_value == pytest.approx(mean)
    assert d2.std_dev == pytest.approx(std_dev)
    assert d2.lower == pytest.approx(lower)
    assert d2.upper == pytest.approx(upper)

    # Test PDF evaluation
    # PDF should be zero outside bounds
    assert d.evaluate(lower - 1.0) == 0.0
    assert d.evaluate(upper + 1.0) == 0.0

    # PDF should be positive inside bounds
    assert d.evaluate(mean) > 0.0

    # PDF should be higher than untruncated at the mean (due to renormalization)
    d_unbounded = openmc.stats.Normal(mean, std_dev)
    assert d.evaluate(mean) > d_unbounded.evaluate(mean)

    # Verify that PDF integrates to approximately 1
    x = np.linspace(lower, upper, 1000)
    integral = trapezoid(d.evaluate(x), x)
    assert integral == pytest.approx(1.0, rel=0.01)

    # Sample truncated distribution
    n_samples = 10_000
    samples, weights = d.sample(n_samples)

    # All samples should be within bounds
    assert np.all(samples >= lower)
    assert np.all(samples <= upper)

    # Weights should all be 1 (no biasing)
    assert np.all(weights == 1.0)


def test_normal_truncated_one_sided():
    # Test lower-bounded only (positive half-normal centered at 0)
    d_lower = openmc.stats.Normal(0.0, 1.0, lower=0.0)
    assert d_lower.lower == 0.0
    assert d_lower.upper == np.inf
    assert d_lower.evaluate(-1.0) == 0.0
    assert d_lower.evaluate(1.0) > 0.0

    # PDF at 0 should be approximately 2 * 0.3989 ≈ 0.798 (half-normal)
    assert d_lower.evaluate(0.0) == pytest.approx(0.798, rel=0.01)

    # Test upper-bounded only
    d_upper = openmc.stats.Normal(0.0, 1.0, upper=0.0)
    assert d_upper.lower == -np.inf
    assert d_upper.upper == 0.0
    assert d_upper.evaluate(1.0) == 0.0
    assert d_upper.evaluate(-1.0) > 0.0


def test_normal_truncated_errors():
    # Invalid bounds (lower >= upper)
    with pytest.raises(ValueError):
        openmc.stats.Normal(0.0, 1.0, lower=1.0, upper=0.0)

    with pytest.raises(ValueError):
        openmc.stats.Normal(0.0, 1.0, lower=1.0, upper=1.0)


@pytest.mark.flaky(reruns=1)
def test_muir():
    mean = 10.0
    mass = 5.0
    temp = 20000.
    d = openmc.stats.muir(mean, mass, temp)
    assert isinstance(d, openmc.stats.Normal)

    elem = d.to_xml_element('energy')
    assert elem.attrib['type'] == 'normal'

    d = openmc.stats.Univariate.from_xml_element(elem)
    assert isinstance(d, openmc.stats.Normal)

    # sample muir distribution
    n_samples = 100_000
    samples, weights = d.sample(n_samples)
    assert_sample_mean(samples, mean)
    assert np.all(weights == 1.0)


@pytest.mark.flaky(reruns=1)
def test_combine_distributions():
    # Combine two discrete (same data as in test_merge_discrete)
    x1 = [0.0, 1.0, 10.0]
    p1 = [0.3, 0.2, 0.5]
    d1 = openmc.stats.Discrete(x1, p1)
    x2 = [0.5, 1.0, 5.0]
    p2 = [0.4, 0.5, 0.1]
    d2 = openmc.stats.Discrete(x2, p2)

    # Merged distribution should have x values sorted and probabilities
    # appropriately combined. Duplicate x values should appear once.
    merged = openmc.stats.combine_distributions([d1, d2], [0.6, 0.4])
    assert isinstance(merged, openmc.stats.Discrete)
    assert merged.x == pytest.approx([0.0, 0.5, 1.0, 5.0, 10.0])
    assert merged.p == pytest.approx(
        [0.6*0.3, 0.4*0.4, 0.6*0.2 + 0.4*0.5, 0.4*0.1, 0.6*0.5])

    # Probabilities add up but are not normalized
    d1 = openmc.stats.Discrete([3.0], [1.0])
    triple = openmc.stats.combine_distributions([d1, d1, d1], [1.0, 2.0, 3.0])
    assert triple.x == pytest.approx([3.0])
    assert triple.p == pytest.approx([6.0])

    # Combine discrete and tabular
    t1 = openmc.stats.Tabular(x2, p2)
    mixed = openmc.stats.combine_distributions([d1, t1], [0.5, 0.5])
    assert isinstance(mixed, openmc.stats.Mixture)
    assert len(mixed.distribution) == 2
    assert len(mixed.probability) == 2
    assert mixed == openmc.stats.combine_distributions([mixed], [1.0])

    # Mixture combined with another distribution: probabilities should be
    # correctly scaled when the Mixture is flattened
    d_a = openmc.stats.delta_function(1.0)
    d_b = openmc.stats.delta_function(2.0)
    m = openmc.stats.Mixture([0.3, 0.7], [d_a, d_b])
    extra = openmc.stats.delta_function(3.0)
    result = openmc.stats.combine_distributions([m, extra], [0.5, 0.5])
    assert isinstance(result, openmc.stats.Discrete)
    assert result.x == pytest.approx([1.0, 2.0, 3.0])
    assert result.p == pytest.approx([0.5*0.3, 0.5*0.7, 0.5])

    # Passing a Mixture with a bias should warn that the bias is dropped
    biased_m = openmc.stats.Mixture([0.5, 0.5], [d_a, d_b], bias=[0.8, 0.2])
    with pytest.warns(UserWarning, match='bias'):
        openmc.stats.combine_distributions([biased_m], [1.0])

    # Single tabular returns a tabular distribution with scaled probabilities
    t_single = openmc.stats.Tabular([0.0, 1.0], [2.0, 0.0])
    scaled = openmc.stats.combine_distributions([t_single], [0.25])
    assert isinstance(scaled, openmc.stats.Tabular)
    assert scaled.p == pytest.approx([0.5, 0.0])

    # Mixture with biased tabular should preserve unbiased mean via weights
    bias = openmc.stats.Tabular([0.0, 1.0], [2.0, 0.0])
    t_biased = openmc.stats.Tabular([0.0, 1.0], [1.0, 1.0], bias=bias)
    d1 = openmc.stats.delta_function(0.0)
    mixed = openmc.stats.combine_distributions([t_biased, d1], [0.5, 0.5])
    assert isinstance(mixed, openmc.stats.Mixture)
    samples, weights = mixed.sample(10_000)
    assert_sample_mean(samples*weights, 0.25)

    # Combine 1 discrete and 2 tabular -- the tabular distributions should
    # combine to produce a uniform distribution with mean 0.5. The combined
    # distribution should have a mean of 0.25.
    t1 = openmc.stats.Tabular([0., 1.], [2.0, 0.0])
    t2 = openmc.stats.Tabular([0., 1.], [0.0, 2.0])
    d1 = openmc.stats.delta_function(0.0)
    combined = openmc.stats.combine_distributions([t1, t2, d1], [0.25, 0.25, 0.5])
    assert combined.integral() == pytest.approx(1.0)

    # Sample the combined distribution and make sure the sample mean is within
    # uncertainty of the expected value
    samples, weights = combined.sample(10_000)
    assert_sample_mean(samples, 0.25)
    assert np.all(weights == 1.0)

    # If biased/unbiased Discrete distributions are combined, unbiased probability
    # should be conserved and points from both original distributions should be
    # assigned bias probabilities.
    x1 = [0.0, 1.0, 10.0]
    p1 = [0.3, 0.2, 0.5]
    b1 = [0.2, 0.5, 0.3]
    d1 = openmc.stats.Discrete(x1, p1, b1)
    x2 = [0.5, 1.0, 5.0]
    p2 = [0.4, 0.5, 0.1]
    d2 = openmc.stats.Discrete(x2, p2)
    combined = openmc.stats.combine_distributions([d1, d2, t1], [0.25, 0.25, 0.5])

    p3 = [0.075, 0.1, 0.175, 0.025, 0.125]
    b3 = [0.05, 0.1, 0.25, 0.025, 0.075]
    assert all(combined.distribution[-1].p == p3)
    assert all(combined.distribution[-1].bias == b3)


def test_reference_vwu_projection():
    """When a non-orthogonal vector is provided, the setter should project out
    any component along reference_uvw so the stored vector is orthogonal.
    """
    pa = openmc.stats.PolarAzimuthal()  # default reference_uvw == (0, 0, 1)

    # Provide a vector that is not orthogonal to (0,0,1)
    pa.reference_vwu = (2.0, 0.5, 0.3)

    reference_v = np.asarray(pa.reference_vwu)
    reference_u = np.asarray(pa.reference_uvw)

    # reference_v should be orthogonal to reference_u
    assert abs(np.dot(reference_v, reference_u)) < 1e-6


def test_reference_vwu_normalization():
    """When a non-normalized vector is provided, the setter should normalize
    the projected vector to unit length.
    """
    pa = openmc.stats.PolarAzimuthal()  # default reference_uvw == (0, 0, 1)

    # Provide a vector that is neither orthogonal to (0,0,1) nor unit-length
    pa.reference_vwu = (2.0, 0.5, 0.3)

    reference_v = np.asarray(pa.reference_vwu)

    # reference_v should be unit length
    assert np.isclose(np.linalg.norm(reference_v), 1.0, atol=1e-12)


def test_fusion_spectrum_dd():
    d = openmc.stats.fusion_neutron_spectrum(10e3, 'DD')
    assert isinstance(d, openmc.stats.Normal)

    # E_0 for D(d,n)3He is ~2.45 MeV; thermal shift at 10 keV should be
    # several tens of keV, so mean should be noticeably above E_0
    assert d.mean_value > 2.45e6
    assert d.mean_value < 2.6e6

    # Standard deviation should be positive and on order of ~50-100 keV
    assert d.std_dev > 30e3
    assert d.std_dev < 200e3


def test_fusion_spectrum_dt():
    d = openmc.stats.fusion_neutron_spectrum(10e3, 'DT')
    assert isinstance(d, openmc.stats.Normal)

    # E_0 for T(d,n)alpha is ~14.02 MeV; with thermal shift mean should be
    # above E_0 by several tens of keV
    assert d.mean_value > 14.02e6
    assert d.mean_value < 14.2e6

    # Standard deviation should be on order of ~200-400 keV
    assert d.std_dev > 100e3
    assert d.std_dev < 500e3


def test_fusion_spectrum_temp_continuity():
    # Verify the low-T and high-T formulas produce nearly identical results
    # at the 40 keV switchover point
    d_lo = openmc.stats.fusion_neutron_spectrum(39.99e3, 'DT')
    d_hi = openmc.stats.fusion_neutron_spectrum(40.01e3, 'DT')

    assert d_lo.mean_value == pytest.approx(d_hi.mean_value, rel=1e-3)
    assert d_lo.std_dev == pytest.approx(d_hi.std_dev, rel=1e-3)

    # Same check for DD
    d_lo = openmc.stats.fusion_neutron_spectrum(39.99e3, 'DD')
    d_hi = openmc.stats.fusion_neutron_spectrum(40.01e3, 'DD')

    assert d_lo.mean_value == pytest.approx(d_hi.mean_value, rel=1e-3)
    assert d_lo.std_dev == pytest.approx(d_hi.std_dev, rel=1e-3)


def test_fusion_spectrum_high_temp():
    # At T_i = 80 keV (high-T regime), ensure the function still produces
    # reasonable results using Table IV formulas
    for reactants in ('DD', 'DT'):
        d = openmc.stats.fusion_neutron_spectrum(80e3, reactants)
        assert isinstance(d, openmc.stats.Normal)
        assert d.mean_value > 0
        assert d.std_dev > 0

    # DT mean at 80 keV should be higher than at 10 keV
    d_10 = openmc.stats.fusion_neutron_spectrum(10e3, 'DT')
    d_80 = openmc.stats.fusion_neutron_spectrum(80e3, 'DT')
    assert d_80.mean_value > d_10.mean_value
    assert d_80.std_dev > d_10.std_dev


def test_fusion_spectrum_zero_temp():
    # At very low temperature, mean should approach E_0 and width should
    # approach zero
    d = openmc.stats.fusion_neutron_spectrum(1.0, 'DT')
    assert d.mean_value == pytest.approx(14.049e6, rel=1e-3)
    assert d.std_dev < 5e3  # width approaches zero at low temperature


def test_fusion_spectrum_invalid():
    # Invalid reactant string should raise an error
    with pytest.raises(ValueError):
        openmc.stats.fusion_neutron_spectrum(10e3, '🐔🧇')

    # Negative temperature should raise an error
    with pytest.raises(ValueError):
        openmc.stats.fusion_neutron_spectrum(-10e3, 'DT')

    # Temperature above 100 keV should raise an error
    with pytest.raises(ValueError):
        openmc.stats.fusion_neutron_spectrum(101e3, 'DT')


@pytest.fixture(autouse=False)
def decay_spectrum_chain():
    """Set chain_file for the duration of a test and clear the _photon_integral
    cache so results from a different chain don't bleed across tests."""
    CHAIN_FILE = (Path(__file__).parents[1] / 'chain_simple.xml').resolve()
    DecaySpectrum._photon_integral.cache_clear()
    with openmc.config.patch('chain_file', CHAIN_FILE):
        yield
    DecaySpectrum._photon_integral.cache_clear()


def test_decay_spectrum_construction():
    nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4}
    d = openmc.stats.DecaySpectrum(nuclides, volume=100.0)
    assert d.nuclides == nuclides
    assert d.volume == pytest.approx(100.0)
    assert len(d) == 2


def test_decay_spectrum_validation():
    # nuclides must be a dict
    with pytest.raises(TypeError):
        openmc.stats.DecaySpectrum(['I135'], volume=1.0)

    # densities must be > 0
    with pytest.raises(ValueError):
        openmc.stats.DecaySpectrum({'I135': -1.0}, volume=1.0)

    # volume must be > 0
    with pytest.raises(ValueError):
        openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=-1.0)

    with pytest.raises(ValueError):
        openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=0.0)


def test_decay_spectrum_xml_roundtrip():
    nuclides = {'I135': 1.5e-3, 'Xe135': 8.2e-4}
    d = openmc.stats.DecaySpectrum(nuclides, volume=100.0)

    elem = d.to_xml_element('energy')
    assert elem.get('type') == 'decay_spectrum'
    assert float(elem.get('volume')) == pytest.approx(100.0)
    assert elem.findtext('nuclides').split() == list(nuclides)
    assert [float(x) for x in elem.findtext('parameters').split()] == pytest.approx(
        list(nuclides.values()))

    # Round-trip via DecaySpectrum.from_xml_element
    d2 = openmc.stats.DecaySpectrum.from_xml_element(elem)
    assert d2.nuclides == nuclides
    assert d2.volume == pytest.approx(100.0)

    # Round-trip via the Univariate dispatcher
    d3 = openmc.stats.Univariate.from_xml_element(elem)
    assert isinstance(d3, openmc.stats.DecaySpectrum)
    assert d3 == d


def test_decay_spectrum_to_distribution(decay_spectrum_chain):
    # Single emitting nuclide -> concrete distribution, not None
    d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0)
    dist = d.to_distribution()
    assert dist is not None

    # Result is cached on second call
    dist2 = d.to_distribution()
    assert dist2 is dist

    # Nuclide with no photon source -> None
    d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0)
    assert d_stable.to_distribution() is None

    # Mixture of emitters -> non-None combined distribution
    d_mix = openmc.stats.DecaySpectrum(
        {'I135': 1e-3, 'Xe135': 5e-4}, volume=10.0
    )
    dist_mix = d_mix.to_distribution()
    assert dist_mix is not None


def test_decay_spectrum_integral(decay_spectrum_chain):
    # For an emitting nuclide, integral should be > 0
    d = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=10.0)
    assert d.integral() > 0.0

    # Proportional to density: doubling density doubles integral
    d2 = openmc.stats.DecaySpectrum({'I135': 2e-3}, volume=10.0)
    assert d2.integral() == pytest.approx(2.0 * d.integral())

    # Proportional to volume
    d3 = openmc.stats.DecaySpectrum({'I135': 1e-3}, volume=20.0)
    assert d3.integral() == pytest.approx(2.0 * d.integral())

    # Pure non-emitter -> 0.0
    d_stable = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0)
    assert d_stable.integral() == pytest.approx(0.0)


def test_decay_spectrum_clip(decay_spectrum_chain):
    # Stable / non-emitting nuclides are removed unconditionally
    d = openmc.stats.DecaySpectrum(
        {'I135': 1e-3, 'Xe135': 5e-4, 'Xe136': 1.0, 'Cs135': 1.0},
        volume=10.0,
    )
    d_clip = d.clip()
    assert 'Xe136' not in d_clip.nuclides
    assert 'Cs135' not in d_clip.nuclides
    assert 'I135' in d_clip.nuclides
    assert 'Xe135' in d_clip.nuclides
    # Original is unchanged
    assert 'Xe136' in d.nuclides

    # inplace=True modifies and returns the same object
    d_same = d.clip(inplace=True)
    assert d_same is d
    assert 'Xe136' not in d.nuclides

    # A nuclide with negligible emission rate is removed by tolerance clipping.
    # U235 has a very small integral (~4e-17 Bq/atom) compared with I135 (~4e-5)
    d_tight = openmc.stats.DecaySpectrum(
        {'I135': 1e-3, 'U235': 1e-3}, volume=10.0
    )
    d_tight_clip = d_tight.clip(tolerance=1e-9)
    assert 'U235' not in d_tight_clip.nuclides
    assert 'I135' in d_tight_clip.nuclides

    # All non-emitters -> empty nuclides dict
    d_empty = openmc.stats.DecaySpectrum({'Xe136': 1e-3}, volume=10.0)
    d_empty.clip(inplace=True)
    assert d_empty.nuclides == {}