OpenMC/openmc/weight_windows.py

from __future__ import annotations
from numbers import Real, Integral
from collections.abc import Iterable, Sequence
from pathlib import Path
from typing import Self
import warnings

import lxml.etree as ET
import numpy as np
import h5py

import openmc
from openmc.mesh import MeshBase, RectilinearMesh, CylindricalMesh, SphericalMesh, UnstructuredMesh
from openmc.tallies import Tallies
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from ._xml import get_elem_list, get_text, clean_indentation
from .mixin import IDManagerMixin
from .particle_type import ParticleType


class WeightWindows(IDManagerMixin):
    """Mesh-based weight windows

    This class enables you to specify weight window parameters that are used in
    a simulation. Multiple sets of weight windows can be defined for different
    meshes and different particles. An iterable of :class:`WeightWindows`
    instances can be assigned to the :attr:`openmc.Settings.weight_windows`
    attribute, which is then exported to XML.

    Weight window lower/upper bounds are to be specified for each combination of
    a mesh element and an energy bin. Thus the total number of bounds should be
    equal to the product of the number of mesh bins and the number of energy
    bins.

    .. versionadded:: 0.13

    Parameters
    ----------
    mesh : openmc.MeshBase
        Mesh for the weight windows
    lower_ww_bounds : Iterable of Real
        A list of values for which each value is the lower bound of a weight
        window
    upper_ww_bounds : Iterable of Real
        A list of values for which each value is the upper bound of a weight
        window
    upper_bound_ratio : float
        Ratio of the lower to upper weight window bounds
    energy_bounds : Iterable of Real
        A list of values for which each successive pair constitutes a range of
        energies in [eV] for a single bin. If no energy bins are provided, the
        maximum and minimum energy for the data available at runtime.
    particle_type : str or int or openmc.ParticleType
        Particle type the weight windows apply to
    survival_ratio : float
        Ratio of the survival weight to the lower weight window bound for
        rouletting
    max_lower_bound_ratio : float
        Maximum allowed ratio of a particle's weight to the weight window's
        lower bound. A factor will be applied to raise the weight window to be
        lower than the particle's weight by a factor of max_lower_bound_ratio
        during transport if exceeded.
    max_split : int
        Maximum allowable number of particles when splitting
    weight_cutoff : float
        Threshold below which particles will be terminated
    id : int
       Unique identifier for the weight window settings. If not specified, an
       identifier will automatically be assigned.

    Attributes
    ----------
    id : int
       Unique identifier for the weight window settings.
    mesh : openmc.MeshBase
        Mesh for the weight windows with dimension (ni, nj, nk)
    particle_type : str
        Particle type the weight windows apply to
    energy_bounds : Iterable of Real
        A list of values for which each successive pair constitutes a range of
        energies in [eV] for a single bin
    num_energy_bins : int
        Number of energy bins
    lower_ww_bounds : numpy.ndarray of float
        An array of values for which each value is the lower bound of a weight
        window. Shape: (ni, nj, nk, num_energy_bins) for StructuredMesh;
        (num_elements, num_energy_bins) for UnstructuredMesh
    upper_ww_bounds : numpy.ndarray of float
        An array of values for which each value is the upper bound of a weight
        window. Shape: (ni, nj, nk, num_energy_bins) for StructuredMesh;
        (num_elements, num_energy_bins) for UnstructuredMesh
    survival_ratio : float
        Ratio of the survival weight to the lower weight window bound for
        rouletting
    max_lower_bound_ratio : float
        Maximum allowed ratio of a particle's weight to the weight window's
        lower bound. (Default: 1.0)
    max_split : int
        Maximum allowable number of particles when splitting
    weight_cutoff : float
        Threshold below which particles will be terminated

    See Also
    --------
    openmc.Settings

    """
    next_id = 1
    used_ids = set()

    def __init__(
        self,
        mesh: MeshBase,
        lower_ww_bounds: Iterable[float],
        upper_ww_bounds: Iterable[float] | None = None,
        upper_bound_ratio: float | None = None,
        energy_bounds: Iterable[Real] | None = None,
        particle_type: str | int | openmc.ParticleType = 'neutron',
        survival_ratio: float = 3.0,
        max_lower_bound_ratio: float | None = None,
        max_split: int = 10,
        weight_cutoff: float = 1.e-38,
        id: int | None = None
    ):
        self.mesh = mesh
        self.id = id
        self.particle_type = particle_type
        self._energy_bounds = None
        if energy_bounds is not None:
            self.energy_bounds = energy_bounds
        self.lower_ww_bounds = lower_ww_bounds

        if upper_ww_bounds is not None and upper_bound_ratio:
            raise ValueError("Exactly one of upper_ww_bounds and "
                             "upper_bound_ratio must be present.")

        if upper_ww_bounds is None and upper_bound_ratio is None:
            raise ValueError("Exactly one of upper_ww_bounds and "
                             "upper_bound_ratio must be present.")

        if upper_bound_ratio:
            self.upper_ww_bounds = [
                lb * upper_bound_ratio for lb in self.lower_ww_bounds
            ]

        if upper_ww_bounds is not None:
            self.upper_ww_bounds = upper_ww_bounds

        if len(self.lower_ww_bounds) != len(self.upper_ww_bounds):
            raise ValueError('Size of the lower and upper weight '
                             'window bounds do not match')

        self.survival_ratio = survival_ratio

        self._max_lower_bound_ratio = None
        if max_lower_bound_ratio is not None:
            self.max_lower_bound_ratio = max_lower_bound_ratio

        self.max_split = max_split
        self.weight_cutoff = weight_cutoff

    def __repr__(self) -> str:
        string = type(self).__name__ + '\n'
        string += '{: <16}=\t{}\n'.format('\tID', self._id)
        string += '{: <16}=\t{}\n'.format('\tMesh', self.mesh)
        string += '{: <16}=\t{}\n'.format('\tParticle Type', self._particle_type)
        string += '{: <16}=\t{}\n'.format('\tEnergy Bounds', self._energy_bounds)
        string += '{: <16}=\t{}\n'.format('\tMax lower bound ratio', self.max_lower_bound_ratio)
        string += '{: <16}=\t{}\n'.format('\tLower WW Bounds', self._lower_ww_bounds)
        string += '{: <16}=\t{}\n'.format('\tUpper WW Bounds', self._upper_ww_bounds)
        string += '{: <16}=\t{}\n'.format('\tSurvival Ratio', self._survival_ratio)
        string += '{: <16}=\t{}\n'.format('\tMax Split', self._max_split)
        string += '{: <16}=\t{}\n'.format('\tWeight Cutoff', self._weight_cutoff)
        return string

    def __eq__(self, other: WeightWindows) -> bool:
        # ensure that `other` is a WeightWindows object
        if not isinstance(other, WeightWindows):
            return False

        # TODO: add ability to check mesh equality

        # check several attributes directly
        attrs = ('particle_type',
                 'survival_ratio',
                 'max_lower_bound_ratio',
                 'max_split',
                 'weight_cutoff')
        for attr in attrs:
            if getattr(self, attr) != getattr(other, attr):
                return False

        # save most expensive checks for last
        if not np.array_equal(self.energy_bounds, other.energy_bounds):
            return False

        if not np.array_equal(self.lower_ww_bounds, other.lower_ww_bounds):
            return False

        if not np.array_equal(self.upper_ww_bounds, other.upper_ww_bounds):
            return False

        return True

    @property
    def mesh(self) -> MeshBase:
        return self._mesh

    @mesh.setter
    def mesh(self, mesh: MeshBase):
        cv.check_type('Weight window mesh', mesh, MeshBase)
        self._mesh = mesh

    @property
    def particle_type(self) -> ParticleType:
        return self._particle_type

    @particle_type.setter
    def particle_type(self, pt):
        ptype = ParticleType(pt)
        if ptype not in {ParticleType.NEUTRON, ParticleType.PHOTON}:
            raise ValueError("Weight windows can only be applied for neutrons or photons")
        self._particle_type = ptype

    @property
    def energy_bounds(self) -> Iterable[Real]:
        return self._energy_bounds

    @energy_bounds.setter
    def energy_bounds(self, bounds: Iterable[float]):
        cv.check_type('Energy bounds', bounds, Iterable, Real)
        self._energy_bounds = np.asarray(bounds)

    @property
    def num_energy_bins(self) -> int:
        if self.energy_bounds is None:
            return 1
        return self.energy_bounds.size - 1

    @property
    def lower_ww_bounds(self) -> np.ndarray:
        return self._lower_ww_bounds

    @lower_ww_bounds.setter
    def lower_ww_bounds(self, bounds: Iterable[float]):
        cv.check_iterable_type('Lower WW bounds',
                               bounds,
                               Real,
                               min_depth=1,
                               max_depth=4)
        # reshape data according to mesh and energy bins
        bounds = np.asarray(bounds)
        if isinstance(self.mesh, UnstructuredMesh):
            bounds = bounds.reshape(-1, self.num_energy_bins)
        else:
            bounds = bounds.reshape(*self.mesh.dimension, self.num_energy_bins)
        self._lower_ww_bounds = bounds

    @property
    def upper_ww_bounds(self) -> np.ndarray:
        return self._upper_ww_bounds

    @upper_ww_bounds.setter
    def upper_ww_bounds(self, bounds: Iterable[float]):
        cv.check_iterable_type('Upper WW bounds',
                               bounds,
                               Real,
                               min_depth=1,
                               max_depth=4)
        # reshape data according to mesh and energy bins
        bounds = np.asarray(bounds)
        if isinstance(self.mesh, UnstructuredMesh):
            bounds = bounds.reshape(-1, self.num_energy_bins)
        else:
            bounds = bounds.reshape(*self.mesh.dimension, self.num_energy_bins)
        self._upper_ww_bounds = bounds

    @property
    def survival_ratio(self) -> float:
        return self._survival_ratio

    @survival_ratio.setter
    def survival_ratio(self, val: float):
        cv.check_type('Survival ratio', val, Real)
        cv.check_greater_than('Survival ratio', val, 1.0, True)
        self._survival_ratio = val

    @property
    def max_lower_bound_ratio(self) -> float:
        return self._max_lower_bound_ratio

    @max_lower_bound_ratio.setter
    def max_lower_bound_ratio(self, val: float):
        cv.check_type('Maximum lower bound ratio', val, Real)
        cv.check_greater_than('Maximum lower bound ratio', val, 1.0, equality=True)
        self._max_lower_bound_ratio = val

    @property
    def max_split(self) -> int:
        return self._max_split

    @max_split.setter
    def max_split(self, val: int):
        cv.check_type('Max split', val, Integral)
        self._max_split = val

    @property
    def weight_cutoff(self) -> float:
        return self._weight_cutoff

    @weight_cutoff.setter
    def weight_cutoff(self, cutoff: float):
        cv.check_type('Weight cutoff', cutoff, Real)
        cv.check_greater_than('Weight cutoff', cutoff, 0.0, True)
        self._weight_cutoff = cutoff

    def to_xml_element(self) -> ET.Element:
        """Return an XML representation of the weight window settings

        Returns
        -------
        element : lxml.etree._Element
            XML element containing the weight window information
        """
        element = ET.Element('weight_windows')

        element.set('id', str(self._id))

        subelement = ET.SubElement(element, 'mesh')
        subelement.text = str(self.mesh.id)

        subelement = ET.SubElement(element, 'particle_type')
        subelement.text = str(self.particle_type)

        if self.energy_bounds is not None:
            subelement = ET.SubElement(element, 'energy_bounds')
            subelement.text = ' '.join(str(e) for e in self.energy_bounds)

        subelement = ET.SubElement(element, 'lower_ww_bounds')
        subelement.text = ' '.join(str(b) for b in self.lower_ww_bounds.ravel('F'))

        subelement = ET.SubElement(element, 'upper_ww_bounds')
        subelement.text = ' '.join(str(b) for b in self.upper_ww_bounds.ravel('F'))

        subelement = ET.SubElement(element, 'survival_ratio')
        subelement.text = str(self.survival_ratio)

        if self.max_lower_bound_ratio is not None:
            subelement = ET.SubElement(element, 'max_lower_bound_ratio')
            subelement.text = str(self.max_lower_bound_ratio)

        subelement = ET.SubElement(element, 'max_split')
        subelement.text = str(self.max_split)

        subelement = ET.SubElement(element, 'weight_cutoff')
        subelement.text = str(self.weight_cutoff)

        return element

    @classmethod
    def from_xml_element(cls, elem: ET.Element, meshes: dict[int, MeshBase]) -> Self:
        """Generate weight window settings from an XML element

        Parameters
        ----------
        elem : lxml.etree._Element
            XML element
        meshes : dict
            Dictionary mapping IDs to mesh objects

        Returns
        -------
        openmc.WeightWindows
            Weight windows object
        """
        # Get mesh for weight windows
        mesh_id = int(get_text(elem, 'mesh'))
        if mesh_id not in meshes:
            raise ValueError(f'Could not locate mesh with ID "{mesh_id}"')
        mesh = meshes[mesh_id]

        # Read all other parameters
        lower_ww_bounds = get_elem_list(elem, "lower_ww_bounds", float)
        upper_ww_bounds = get_elem_list(elem, "upper_ww_bounds", float)
        e_bounds = get_elem_list(elem, "energy_bounds", float)
        particle_type = get_text(elem, 'particle_type')
        survival_ratio = float(get_text(elem, 'survival_ratio'))

        ww_shape = (len(e_bounds) - 1,) + mesh.dimension[::-1]
        lower_ww_bounds = np.array(lower_ww_bounds).reshape(ww_shape).T
        upper_ww_bounds = np.array(upper_ww_bounds).reshape(ww_shape).T

        max_lower_bound_ratio = None
        if get_text(elem, 'max_lower_bound_ratio'):
            max_lower_bound_ratio = float(get_text(elem, 'max_lower_bound_ratio'))

        max_split = int(get_text(elem, 'max_split'))
        weight_cutoff = float(get_text(elem, 'weight_cutoff'))
        id = int(get_text(elem, 'id'))

        return cls(
            mesh=mesh,
            lower_ww_bounds=lower_ww_bounds,
            upper_ww_bounds=upper_ww_bounds,
            energy_bounds=e_bounds,
            particle_type=particle_type,
            survival_ratio=survival_ratio,
            max_lower_bound_ratio=max_lower_bound_ratio,
            max_split=max_split,
            weight_cutoff=weight_cutoff,
            id=id
        )

    @classmethod
    def from_hdf5(cls, group: h5py.Group, meshes: dict[int, MeshBase]) -> Self:
        """Create weight windows from HDF5 group

        Parameters
        ----------
        group : h5py.Group
            Group in HDF5 file
        meshes : dict
            Dictionary mapping IDs to mesh objects

        Returns
        -------
        openmc.WeightWindows
            A weight window object
        """

        id = int(group.name.split('/')[-1].lstrip('weight_windows'))
        mesh_id = group['mesh'][()]
        mesh = meshes[mesh_id]

        ptype = group['particle_type'][()].decode()
        e_bounds = group['energy_bounds'][()]
        # weight window bounds are stored with the shape (e, k, j, i)
        # in C++ and HDF5 -- the opposite of how they are stored here
        shape = (e_bounds.size - 1,  *mesh.dimension[::-1])
        lower_ww_bounds = group['lower_ww_bounds'][()].reshape(shape).T
        upper_ww_bounds = group['upper_ww_bounds'][()].reshape(shape).T
        survival_ratio = group['survival_ratio'][()]

        max_lower_bound_ratio = None
        if group.get('max_lower_bound_ratio') is not None:
            max_lower_bound_ratio = group['max_lower_bound_ratio'][()]

        max_split = group['max_split'][()]
        weight_cutoff = group['weight_cutoff'][()]

        return cls(
            mesh=mesh,
            lower_ww_bounds=lower_ww_bounds,
            upper_ww_bounds=upper_ww_bounds,
            energy_bounds=e_bounds,
            particle_type=ptype,
            survival_ratio=survival_ratio,
            max_lower_bound_ratio=max_lower_bound_ratio,
            max_split=max_split,
            weight_cutoff=weight_cutoff,
            id=id
        )


def wwinp_to_wws(path: PathLike) -> WeightWindowsList:
    """Create WeightWindows instances from a wwinp file

    .. versionadded:: 0.13.1

    Parameters
    ----------
    path : str or pathlib.Path
        Path to the wwinp file

    Returns
    -------
    WeightWindowsList
    """
    warnings.warn(
        "This function is deprecated in favor of 'WeightWindowsList.from_wwinp'",
        FutureWarning
    )
    return WeightWindowsList.from_wwinp(path)


class WeightWindowGenerator:
    """Class passed to setting to govern weight window generation
    using the OpenMC executable

    Parameters
    ----------
    mesh : :class:`openmc.MeshBase`
        Mesh used to represent the weight windows spatially
    energy_bounds : Iterable of Real
        A list of values for which each successive pair constitutes a range of
        energies in [eV] for a single bin. If no energy bins are provided, the
        maximum and minimum energy for the data available at runtime.
    particle_type : str or int or openmc.ParticleType
        Particle type the weight windows apply to
    method : {'magic', 'fw_cadis'}
        The weight window generation methodology applied during an update.
    targets : :class:`openmc.Tallies` or iterable of int
        Target tallies for local variance reduction via FW-CADIS.
    max_realizations : int
        The upper limit for number of tally realizations when generating weight
        windows.
    update_interval : int
        The number of tally realizations between updates.
    on_the_fly : bool
        Whether or not to apply weight windows on the fly.

    Attributes
    ----------
    mesh : openmc.MeshBase
        Mesh used to represent the weight windows spatially
    energy_bounds : Iterable of Real
        A list of values for which each successive pair constitutes a range of
        energies in [eV] for a single bin
    particle_type : openmc.ParticleType
        Particle type the weight windows apply to
    method : {'magic', 'fw_cadis'}
        The weight window generation methodology applied during an update.
    targets : :class:`openmc.Tallies` or numpy.ndarray
        Target tallies for local variance reduction via FW-CADIS.
    max_realizations : int
        The upper limit for number of tally realizations when generating weight
        windows.
    update_interval : int
        The number of tally realizations between updates.
    update_parameters : dict
        A set of parameters related to the update.
    on_the_fly : bool
        Whether or not to apply weight windows on the fly.
    """

    _WWG_PARAMS = {'value': str, 'threshold': float, 'ratio': float}

    def __init__(
        self,
        mesh: openmc.MeshBase,
        energy_bounds: Sequence[float] | None = None,
        particle_type: str | int | openmc.ParticleType = 'neutron',
        method: str = 'magic',
        targets: openmc.Tallies | Iterable[int] | None = None,
        max_realizations: int = 1,
        update_interval: int = 1,
        on_the_fly: bool = True
    ):
        self._update_parameters = None

        self.mesh = mesh
        self._energy_bounds = None
        if energy_bounds is not None:
            self.energy_bounds = energy_bounds
        self.particle_type = particle_type
        self.method = method
        self.targets = targets
        self.max_realizations = max_realizations
        self.update_interval = update_interval
        self.on_the_fly = on_the_fly

    def __repr__(self):
        string = type(self).__name__ + '\n'
        string += f'\t{"Mesh":<20}=\t{self.mesh.id}\n'
        string += f'\t{"Particle:":<20}=\t{str(self.particle_type)}\n'
        string += f'\t{"Energy Bounds:":<20}=\t{self.energy_bounds}\n'
        string += f'\t{"Method":<20}=\t{self.method}\n'
        string += f'\t{"Max Realizations:":<20}=\t{self.max_realizations}\n'
        string += f'\t{"Update Interval:":<20}=\t{self.update_interval}\n'
        string += f'\t{"On The Fly:":<20}=\t{self.on_the_fly}\n'
        if self.update_parameters is not None:
            string += f'\t{"Update Parameters:":<20}\n\t\t\t{self.update_parameters}\n'
        string

        return string

    @property
    def mesh(self) -> openmc.MeshBase:
        return self._mesh

    @mesh.setter
    def mesh(self, m: openmc.MeshBase):
        cv.check_type('mesh', m, openmc.MeshBase)
        self._mesh = m

    @property
    def energy_bounds(self) -> Iterable[Real]:
        return self._energy_bounds

    @energy_bounds.setter
    def energy_bounds(self, eb: Iterable[float]):
        cv.check_type('energy bounds', eb, Iterable, Real)
        self._energy_bounds = eb

    @property
    def particle_type(self) -> ParticleType:
        return self._particle_type

    @particle_type.setter
    def particle_type(self, pt):
        ptype = ParticleType(pt)
        if ptype not in {ParticleType.NEUTRON, ParticleType.PHOTON}:
            raise ValueError("Weight windows can only be applied for neutrons or photons")
        self._particle_type = ptype

    @property
    def method(self) -> str:
        return self._method

    @method.setter
    def method(self, m: str):
        cv.check_type('generation method', m, str)
        cv.check_value('generation method', m, ('magic', 'fw_cadis'))
        self._method = m
        if self._update_parameters is not None:
            try:
                self._check_update_parameters()
            except (TypeError, KeyError):
                warnings.warn(f'Update parameters are invalid for the "{m}" method.')
    
    @property
    def targets(self) -> openmc.Tallies:
        return self._targets

    @targets.setter
    def targets(self, t):
        if t is None:
            self._targets = t
        else:
            cv.check_type('Local FW-CADIS target tallies', t, Iterable)
            cv.check_greater_than('Local FW-CADIS target tallies', len(t), 0)
            if not isinstance(t, openmc.Tallies):
                cv.check_iterable_type('Local FW-CADIS target tallies', t, int)
                t = np.asarray(list(t), dtype=int)
            self._targets = t

    @property
    def max_realizations(self) -> int:
        return self._max_realizations

    @max_realizations.setter
    def max_realizations(self, m: int):
        cv.check_type('max tally realizations', m, Integral)
        cv.check_greater_than('max tally realizations', m, 0)
        self._max_realizations = m

    @property
    def update_interval(self) -> int:
        return self._update_interval

    @update_interval.setter
    def update_interval(self, ui: int):
        cv.check_type('update interval', ui, Integral)
        cv.check_greater_than('update interval', ui , 0)
        self._update_interval = ui

    @property
    def update_parameters(self) -> dict:
        return self._update_parameters

    def _check_update_parameters(self, params: dict):
        if self.method == 'magic' or self.method == 'fw_cadis':
            check_params = self._WWG_PARAMS

        for key, val in params.items():
            if key not in check_params:
                raise ValueError(f'Invalid param "{key}" for {self.method} '
                                  'weight window generation')
            cv.check_type(f'weight window generation param: "{key}"', val, self._WWG_PARAMS[key])

    @update_parameters.setter
    def update_parameters(self, params: dict):
        self._check_update_parameters(params)
        self._update_parameters = params

    @property
    def on_the_fly(self) -> bool:
        return self._on_the_fly

    @on_the_fly.setter
    def on_the_fly(self, otf: bool):
        cv.check_type('on the fly generation', otf, bool)
        self._on_the_fly = otf

    def _update_parameters_subelement(self, element: ET.Element):
        if not self.update_parameters:
            return
        params_element = ET.SubElement(element, 'update_parameters')
        for pname, value in self.update_parameters.items():
            param_element = ET.SubElement(params_element, pname)
            param_element.text = str(value)

    @classmethod
    def _sanitize_update_parameters(cls, method: str, update_parameters: dict):
        """
        Attempt to convert update parameters to their appropriate types

        Parameters
        ----------
        method : str
            The update method for which these update parameters should comply
        update_parameters : dict
            The update parameters as-read from the XML node (keys: str, values: str)
        """
        if method == 'magic' or method == 'fw_cadis':
            check_params = cls._WWG_PARAMS

        for param, param_type in check_params.items():
            if param in update_parameters:
                update_parameters[param] = param_type(update_parameters[param])

    def to_xml_element(self):
        """Creates a 'weight_window_generator' element to be written to an XML file.
        """
        element = ET.Element('weight_windows_generator')

        mesh_elem = ET.SubElement(element, 'mesh')
        mesh_elem.text = str(self.mesh.id)
        if self.energy_bounds is not None:
            subelement = ET.SubElement(element, 'energy_bounds')
            subelement.text = ' '.join(str(e) for e in self.energy_bounds)
        particle_elem = ET.SubElement(element, 'particle_type')
        particle_elem.text = str(self.particle_type)
        realizations_elem = ET.SubElement(element, 'max_realizations')
        realizations_elem.text = str(self.max_realizations)
        update_interval_elem = ET.SubElement(element, 'update_interval')
        update_interval_elem.text = str(self.update_interval)
        otf_elem = ET.SubElement(element, 'on_the_fly')
        otf_elem.text = str(self.on_the_fly).lower()
        method_elem = ET.SubElement(element, 'method')
        method_elem.text = self.method
        if self.targets is not None:
            if self.method != 'fw_cadis':
                raise ValueError(
                    "FW-CADIS update method is required in order to use " \
                    "target tallies for WeightWindowGenerator.")
            elif isinstance(self.targets, openmc.Tallies):
                raise RuntimeError(
                    "FW-CADIS target tallies must be checked to ensure they are " \
                    "present on model.tallies. Use model.export_to_xml() or " \
                    "model.export_to_model_xml() to link FW-CADIS target tallies.")
            else:
                targets_elem = ET.SubElement(element, 'targets')
                targets_elem.text = ' '.join(str(tally_id) for tally_id in self.targets)

        if self.update_parameters is not None:
            self._update_parameters_subelement(element)

        clean_indentation(element)

        return element

    @classmethod
    def from_xml_element(cls, elem: ET.Element, meshes: dict) -> Self:
        """
        Create a weight window generation object from an XML element

        Parameters
        ----------
        elem : xml.etree.ElementTree.Element
            XML element
        meshes : dict
            A dictionary with IDs as keys and openmc.MeshBase instances as values

        Returns
        -------
        openmc.WeightWindowGenerator
        """

        mesh_id = int(get_text(elem, 'mesh'))
        mesh = meshes[mesh_id]
        
        energy_bounds = get_elem_list(elem, "energy_bounds", float)
        particle_type = get_text(elem, 'particle_type')

        wwg = cls(mesh, energy_bounds, particle_type)

        wwg.max_realizations = int(get_text(elem, 'max_realizations'))
        wwg.update_interval = int(get_text(elem, 'update_interval'))
        wwg.on_the_fly = bool(get_text(elem, 'on_the_fly'))
        wwg.method = get_text(elem, 'method')
        targets_elem = elem.find('targets')
        if targets_elem is not None:
            if wwg.method != 'fw_cadis':
                raise ValueError(
                    "FW-CADIS update method is required in order to use " \
                    "target tallies for WeightWindowGenerator.")
            else:
                wwg.targets = get_elem_list(elem, "targets")

        if elem.find('update_parameters') is not None:
            update_parameters = {}
            params_elem = elem.find('update_parameters')
            for entry in params_elem:
                update_parameters[entry.tag] = entry.text

            cls._sanitize_update_parameters(wwg.method, update_parameters)
            wwg.update_parameters = update_parameters

        return wwg

def hdf5_to_wws(path='weight_windows.h5') -> WeightWindowsList:
    """Create a WeightWindowsList from a weight windows HDF5 file

    .. versionadded:: 0.14.0

    Parameters
    ----------
    path : cv.PathLike
        Path to the weight windows hdf5 file

    Returns
    -------
    WeightWindowsList
    """
    warnings.warn(
        "This function is deprecated in favor of 'WeightWindowsList.from_hdf5'",
        FutureWarning
    )
    return WeightWindowsList.from_hdf5(path)


class WeightWindowsList(list):
    """A list of WeightWindows objects.

    .. versionadded:: 0.15.3

    Parameters
    ----------
    iterable : iterable of openmc.WeightWindows
        An iterable of WeightWindows objects to initialize the list with

    """
    def __init__(self, iterable: Iterable[WeightWindows] = ()):
        super().__init__(iterable)

    @classmethod
    def from_hdf5(cls, path: PathLike = 'weight_windows.h5') -> Self:
        """Create WeightWindowsList from a weight windows HDF5 file.

        Parameters
        ----------
        path : PathLike
            Path to the weight windows hdf5 file

        Returns
        -------
        WeightWindowsList
            A list of WeightWindows objects read from the file
        """

        with h5py.File(path) as h5_file:
            # read in all of the meshes in the mesh node
            meshes = {}
            for mesh_group in h5_file['meshes']:
                mesh = MeshBase.from_hdf5(h5_file['meshes'][mesh_group])
                meshes[mesh.id] = mesh
            wws = [
                WeightWindows.from_hdf5(ww, meshes)
                for ww in h5_file['weight_windows'].values()
            ]

        return cls(wws)

    @classmethod
    def from_wwinp(cls, path: PathLike) -> Self:
        """Create WeightWindowsList from a wwinp file.

        Parameters
        ----------
        path : PathLike
            Path to the wwinp file

        Returns
        -------
        WeightWindowsList
            A list of WeightWindows objects read from the file
        """

        with open(path) as wwinp:
            # BLOCK 1
            header = wwinp.readline().split(None, 4)
            # read file type, time-dependence, number of
            # particles, mesh type and problem identifier
            _if, iv, ni, nr = [int(x) for x in header[:4]]

            # header value checks
            if _if != 1:
                raise ValueError(f'Found incorrect file type, if: {_if}')

            if iv > 1:
                # read number of time bins for each particle, 'nt(1...ni)'
                nt = np.fromstring(wwinp.readline(), sep=' ', dtype=int)

                # raise error if time bins are present for now
                raise ValueError('Time-dependent weight windows '
                                'are not yet supported')
            else:
                nt = ni * [1]

            # read number of energy bins for each particle, 'ne(1...ni)'
            ne = np.fromstring(wwinp.readline(), sep=' ', dtype=int)

            # read coarse mesh dimensions and lower left corner
            mesh_description = np.fromstring(wwinp.readline(), sep=' ')
            nfx, nfy, nfz = mesh_description[:3].astype(int)
            xyz0 = mesh_description[3:]

            # read cylindrical and spherical mesh vectors if present
            if nr == 16:
                # read number of coarse bins
                line_arr = np.fromstring(wwinp.readline(), sep=' ')
                ncx, ncy, ncz = line_arr[:3].astype(int)
                # read polar vector (x1, y1, z1)
                xyz1 = line_arr[3:]
                # read azimuthal vector (x2, y2, z2)
                line_arr = np.fromstring(wwinp.readline(), sep=' ')
                xyz2 = line_arr[:3]

                # Get polar and azimuthal axes
                polar_axis = xyz1 - xyz0
                azimuthal_axis = xyz2 - xyz0

                # Check for polar axis other than (0, 0, 1)
                norm = np.linalg.norm(polar_axis)
                if not np.isclose(polar_axis[2]/norm, 1.0):
                    raise NotImplementedError('Polar axis not aligned to z-axis not supported')

                # Check for azimuthal axis other than (1, 0, 0)
                norm = np.linalg.norm(azimuthal_axis)
                if not np.isclose(azimuthal_axis[0]/norm, 1.0):
                    raise NotImplementedError('Azimuthal axis not aligned to x-axis not supported')

                # read geometry type
                nwg = int(line_arr[-1])

            elif nr == 10:
                # read rectilinear data:
                # number of coarse mesh bins and mesh type
                ncx, ncy, ncz, nwg = \
                    np.fromstring(wwinp.readline(), sep=' ').astype(int)
            else:
                raise RuntimeError(f'Invalid mesh description (nr) found: {nr}')

            # read BLOCK 2 and BLOCK 3 data into a single array
            ww_data = np.fromstring(wwinp.read(), sep=' ')

        # extract mesh data from the ww_data array
        start_idx = 0

        # first values in the mesh definition arrays are the first
        # coordinate of the grid
        end_idx = start_idx + 1 + 3 * ncx
        i0, i_vals = ww_data[start_idx], ww_data[start_idx+1:end_idx]
        start_idx = end_idx

        end_idx = start_idx + 1 + 3 * ncy
        j0, j_vals = ww_data[start_idx], ww_data[start_idx+1:end_idx]
        start_idx = end_idx

        end_idx = start_idx + 1 + 3 * ncz
        k0, k_vals = ww_data[start_idx], ww_data[start_idx+1:end_idx]
        start_idx = end_idx

        # mesh consistency checks
        if nr == 16 and nwg == 1 or nr == 10 and nwg != 1:
            raise ValueError(f'Mesh description in header ({nr}) '
                            f'does not match the mesh type ({nwg})')

        if nr == 10 and (xyz0 != (i0, j0, k0)).any():
            raise ValueError(f'Mesh origin in the header ({xyz0}) '
                            f' does not match the origin in the mesh '
                            f' description ({i0, j0, k0})')

        # create openmc mesh object
        grids = []
        mesh_definition = [(i0, i_vals, nfx), (j0, j_vals, nfy), (k0, k_vals, nfz)]
        for grid0, grid_vals, n_pnts in mesh_definition:
            # file spec checks for the mesh definition
            if (grid_vals[2::3] != 1.0).any():
                raise ValueError('One or more mesh ratio value, qx, '
                                'is not equal to one')

            s = int(grid_vals[::3].sum())
            if s != n_pnts:
                raise ValueError(f'Sum of the fine bin entries, {s}, does '
                                f'not match the number of fine bins, {n_pnts}')

            # extend the grid based on the next coarse bin endpoint, px
            # and the number of fine bins in the coarse bin, sx
            intervals = grid_vals.reshape(-1, 3)
            coords = [grid0]
            for sx, px, qx in intervals:
                coords += np.linspace(coords[-1], px, int(sx + 1)).tolist()[1:]

            grids.append(np.array(coords))

        if nwg == 1:
            mesh = RectilinearMesh()
            mesh.x_grid, mesh.y_grid, mesh.z_grid = grids
        elif nwg == 2:
            mesh = CylindricalMesh(
                r_grid=grids[0],
                z_grid=grids[1],
                phi_grid=grids[2],
                origin = xyz0,
            )
        elif nwg == 3:
            mesh = SphericalMesh(
                r_grid=grids[0],
                theta_grid=grids[1],
                phi_grid=grids[2],
                origin = xyz0
            )

        # extract weight window values from array
        wws = cls()
        for ne_i, nt_i, particle_type in zip(ne, nt, ('neutron', 'photon')):
            # no information to read for this particle if
            # either the energy bins or time bins are empty
            if ne_i == 0 or nt_i == 0:
                continue

            if iv > 1:
                # time bins are parsed but unused for now
                end_idx = start_idx + nt_i
                time_bounds = ww_data[start_idx:end_idx]
                np.insert(time_bounds, (0,), (0.0,))
                start_idx = end_idx

            # read energy boundaries
            end_idx = start_idx + ne_i
            energy_bounds = np.insert(ww_data[start_idx:end_idx], (0,), (0.0,))
            # convert from MeV to eV
            energy_bounds *= 1e6
            start_idx = end_idx

            # read weight window values
            end_idx = start_idx + (nfx * nfy * nfz) * nt_i * ne_i

            # read values and reshape according to ordering
            # slowest to fastest: t, e, z, y, x
            # reorder with transpose since our ordering is x, y, z, e, t
            ww_shape = (nt_i, ne_i, nfz, nfy, nfx)
            ww_values = ww_data[start_idx:end_idx].reshape(ww_shape).T
            # Only use first time bin since we don't support time dependent weight
            # windows yet.
            ww_values = ww_values[:, :, :, :, 0]
            start_idx = end_idx

            # create a weight window object
            ww = WeightWindows(id=None,
                            mesh=mesh,
                            lower_ww_bounds=ww_values,
                            upper_bound_ratio=5.0,
                            energy_bounds=energy_bounds,
                            particle_type=particle_type)
            wws.append(ww)

        return wws

    def export_to_hdf5(self, path: PathLike = 'weight_windows.h5', **init_kwargs):
        """Write weight windows to an HDF5 file.

        Parameters
        ----------
        path : PathLike
            Path to the file to write weight windows to
        **init_kwargs
            Keyword arguments passed to :func:`openmc.lib.init`

        """
        import openmc.lib
        cv.check_type('path', path, PathLike)

        # Create a temporary model with the weight windows
        model = openmc.Model()
        sph = openmc.Sphere(boundary_type='vacuum')
        cell = openmc.Cell(region=-sph)
        model.geometry = openmc.Geometry([cell])
        model.settings.weight_windows = self
        model.settings.particles = 100
        model.settings.batches = 1

        # Get absolute path before moving to temporary directory
        path = Path(path).resolve()

        # Load the model with openmc.lib and then export it to an HDF5 file
        with openmc.lib.TemporarySession(model, **init_kwargs):
            openmc.lib.export_weight_windows(path)