Repository: terrapower/armi
Pull request #1726 (2024-06-11 01:16:15)

Isolating jagged arrays for numpy compatibility
## What is the change?

Implement a `JaggedArray` class to interface between our jagged data structures (e.g., pin-level data within a block) and HDF5 database without using numpy to flatten the jagged array.

## Why is the change being made?

The current code feeds jagged arrays to numpy and then uses numpy to flatten them for storing in the database and reconstitute them when reading data back from the database. Numpy no longer supports jagged arrays, so we are currently pinned to version 1.23.5 (the latest release is 1.26.4).

The current version of numpy we are using gives a `VisibleDeprecationWarning` about this:

```
.\armi\bookkeeping\db\database3.py:1361: VisibleDeprecationWarning: Creating an ndarray from ragged nested
sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If
you meant to do this, you must specify 'dtype=object' when creating the ndarray.
  data = numpy.array(layout.location)[layoutIndicesForType][
```

This change introduces a `JaggedArray` class that handles the packing into and unpacking out of 1D data structures for the HDF5 database. This allows us to use a current version of numpy.

There is a discussion about this topic: https://github.com/terrapower/armi/discussions/1669

Also, I found this open ticket (https://github.com/terrapower/armi/issues/246) _after_ having done all of the work in this PR. I didn't realize how much others had already thought about how to work around this. Well, I guess I'm throwing my potential solution in the ring. As was noted many times in that ticket, matching the performance of numpy might be the main challenge. Hopefully, the amount of jagged data we have is small enough that the performance hit isn't too significant.

Notes about numpy compatibility
-------------------------------------

Using a current version of numpy, if you try to pass jagged data into `numpy.array()`, you get an error: 

```python
>>> a = [1, 2, 3]
>>> b = [4, 5, 6, 7]
>>> c = [0, 0, 1]
>>> import numpy
>>> d = numpy.array((a, c))
>>> e = numpy.array((b, d))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (2,) + inhomogeneous part.
```

You can get numpy to accept it if you pass in `dtype=object`, but this data type doesn't have the features needed to facilitate flattening the array to store in HDF5, i.e., `array.flatten()` doesn't do anything. So I don't know whether this would help us at all. I'm not 100% positive, but I think we'd still be in the same position of needing to implement the packing and unpacking routines ourselves.

```python
>>> e = numpy.array((b, d), dtype=object)
>>> print(e)
[list([4, 5, 6, 7]) array([[1, 2, 3],
                           [0, 0, 1]])]
>>> e.flatten()
array([list([4, 5, 6, 7]), array([[1, 2, 3],
                                  [0, 0, 1]])], dtype=object)
```

---

## Checklist

- [x] This PR has only [one purpose or idea](https://terrapower.github.io/armi/developer/tooling.html#one-idea-one-pr).
- [x] [Tests](https://terrapower.github.io/armi/developer/tooling.html#test-it) have been added/updated to verify any new/changed code.

<!-- Check the code quality -->

- [x] The code style follows [good practices](https://terrapower.github.io/armi/developer/standards_and_practices.html).
- [x] The commit message(s) follow [good practices](https://terrapower.github.io/armi/developer/tooling.html).

<!-- Check the project-level cruft -->

- [x] The [release notes](https://terrapower.github.io/armi/developer/tooling.html#add-release-notes) have been updated if necessary.
- [x] The [documentation](https://terrapower.github.io/armi/developer/tooling.html#document-it) is still up-to-date in the `doc` folder.
- [x] The dependencies are still up-to-date in `pyproject.toml`.

john-science: I love this! Thanks for getting started!

Can we close this until the PR until you are ready for review?  https://terrapower.github.io/armi/developer/tooling.html#don-t-open-until-it-is-ready

Thanks!
john-science: @mgjarrett  I'm really excited about this PR.  Holler if you need anything.
john-science: @mgjarrett I hope you don't mind, I fixed the merge conflicts in this branch, so the CI would run again.

And while I was at it, I moved `JaggedArray` to it's own file. I can't help myself, sorry!
mgjarrett: > @mgjarrett I hope you don't mind, I fixed the merge conflicts in this branch, so the CI would run again.
> 
> And while I was at it, I moved `JaggedArray` to it's own file. I can't help myself, sorry!

Thank you! I'm hoping to get back to this today, look at your comments and changes and clean up whatever needs to be fixed.
mgjarrett: I did some simple profiling to make sure there's not an unreasonable increase in run time between the old numpy version and the new implementation. (A small increase is probably acceptable because this will enable us to move to Python 3.11 or 3.12, which is reportedly significantly faster than Python 3.9).

The test case I created for profiling is quite simple, consisting of two separate runs. For the first run:

1. Load a reactor operator from inputs
2. Populate `b.p.linPowByPin` with random numbers on each block (this is a jagged dataset)
3. Write the reactor to a database

For the second run:

1. Load from the database saved in the previous run

The time to write the database is about  11-12 seconds, compared to an overall run time of about 100 s. Most of the database writing is related to writing the layout. The `_writeParams` portion is less than 5 seconds, and the `packSpecialData` portion of the database write is less than 1 second in both cases. The database write is slightly faster in the base case (current nala main, numpy 1.23.5), than the new implementation.

Profiling results for writing to the database
-----------------------------------------------

  | writeToDB (s) | _writeParams (s) | packSpecialData (s)
-- | -- | -- | --
Nala Main | 11.5 | 4.39 | 0.336
Jagged Arrays (PR#589) | 11.8 | 4.72 | 0.503
Jagged Arrays (PR#589) (numpy 1.23.5) | 11.9 | 4.97 | 0.559

Profiling results for reading from the database
-----------------------------------------------

Again, the new implementation is slightly slower than the old implementation, although not by much (1.4 s vs. 1.5 s). PR#589 is slower across the board, even for parts of the code that are not touched by the `JaggedArray` update. This difference persisted after multiple runs. It is possibly due to the newer version of numpy (1.26.4) is slower than the older version (1.23.5) on Python 3.9, or something else that has changed in the code. We can't compare on Python 3.11 or 3.12 because numpy 1.23.5 is not compatible. Running the new code with numpy 1.23.5 did slightly improve the run time for `unpackSpecialData` and for the `_readParams`.

  | load (s) | _readParams (s) | unpackSpecialData (s)
-- | -- | -- | --
Nala Main Run 1 | 36.3 | 3.78 | 1.45
Nala Main Run 2 | 40.9 | 3.84 | 1.35
Nala Main Run 3 | 41.6 | 4.28 | 1.35
Jagged Arrays (PR#589) Run 1 | 40.9 | 4.04 | 1.53
Jagged Arrays (PR#589) Run 2 | 40.3 | 4.16 | 1.66
Jagged Arrays (PR#589) Run 3 | 41.6 | 3.75 | 1.40
Jagged Arrays (PR#589) (numpy 1.23.5) | 39.2 | 3.62 | 1.31
john-science: @mgjarrett Your timing results are really interesting. I suppose Python 3.9 is the target version now, but 3.11 will be soon.

On the one hand, people might not like the speed decrease, on the other hand we need to be able to unlock modern versions of NumPy and Python. And both of those are blocked by this NumPy issue.
opotowsky: > @mgjarrett Your timing results are really interesting. I suppose Python 3.9 is the target version now, but 3.11 will be soon.
> 
> On the one hand, people might not like the speed decrease, on the other hand we need to be able to unlock modern versions of NumPy and Python. And both of those are blocked by this NumPy issue.

The speed decrease is not even noticeable (considering the normal variation in run times) for an actual case. I'm quite pumped about this!
john-science: @mgjarrett When this PR is ready for review, please move it from Draft to normal. Right now, I can't tell when you think this PR is ready for prime time.

(And again, ARMI discourages opening Draft PRs.)

Patch:
diff --git a/armi/bookkeeping/db/database3.py b/armi/bookkeeping/db/database3.py
index ee67f70d..df34df43 100644
--- a/armi/bookkeeping/db/database3.py
+++ b/armi/bookkeeping/db/database3.py
@@ -43,36 +43,31 @@ import subprocess
 import sys
 from platform import uname
 from typing import (
-    Optional,
-    Tuple,
-    Type,
-    Dict,
     Any,
+    Dict,
+    Generator,
     List,
+    Optional,
     Sequence,
-    Generator,
+    Tuple,
+    Type,
 )
 
 import h5py
 import numpy
 
-from armi import context
-from armi import getApp
-from armi import meta
-from armi import runLog
-from armi import settings
+from armi import context, getApp, meta, runLog, settings
+from armi.bookkeeping.db.jaggedArray import JaggedArray
 from armi.bookkeeping.db.layout import (
-    Layout,
     DB_VERSION,
+    Layout,
     replaceNonesWithNonsense,
     replaceNonsenseWithNones,
 )
-from armi.bookkeeping.db.typedefs import History, Histories
+from armi.bookkeeping.db.typedefs import Histories, History
 from armi.nucDirectory import nuclideBases
 from armi.physics.neutronics.settings import CONF_LOADING_FILE
-from armi.reactor import grids
-from armi.reactor import parameters
-from armi.reactor import systemLayoutInput
+from armi.reactor import grids, parameters, systemLayoutInput
 from armi.reactor.assemblies import Assembly
 from armi.reactor.blocks import Block
 from armi.reactor.components import Component
@@ -652,7 +647,7 @@ class Database3:
         if groupName in self.h5db:
             return self.h5db[groupName]
         else:
-            group = self.h5db.create_group(groupName)
+            group = self.h5db.create_group(groupName, track_order=True)
             group.attrs["cycle"] = r.p.cycle
             group.attrs["timeNode"] = r.p.timeNode
             return group
@@ -868,14 +863,23 @@ class Database3:
 
         return comp
 
-    def _writeParams(self, h5group, comps):
+    def _writeParams(self, h5group, comps) -> tuple:
+        def _getShape(arr: [numpy.ndarray, List, Tuple]):
+            """Get the shape of a numpy.ndarray, list, or tuple."""
+            if isinstance(arr, numpy.ndarray):
+                return arr.shape
+            elif isinstance(arr, (list, tuple)):
+                return (len(arr),)
+            else:
+                return (1,)
+
         c = comps[0]
         groupName = c.__class__.__name__
         if groupName not in h5group:
             # Only create the group if it doesnt already exist. This happens when
             # re-writing params in the same time node (e.g. something changed between
             # EveryNode and EOC)
-            g = h5group.create_group(groupName)
+            g = h5group.create_group(groupName, track_order=True)
         else:
             g = h5group[groupName]
 
@@ -912,41 +916,48 @@ class Database3:
                     attrs[_SERIALIZER_NAME] = paramDef.serializer.__name__
                     attrs[_SERIALIZER_VERSION] = paramDef.serializer.version
                 else:
-                    data = numpy.array(temp)
+                    # check if temp is a jagged array
+                    if any(isinstance(x, (numpy.ndarray, list)) for x in temp):
+                        jagged = len(set([_getShape(x) for x in temp])) != 1
+                    else:
+                        jagged = False
+                    data = (
+                        JaggedArray(temp, paramDef.name)
+                        if jagged
+                        else numpy.array(temp)
+                    )
                     del temp
 
-            # Convert Unicode to byte-string
-            if data.dtype.kind == "U":
-                data = data.astype("S")
-
-            if data.dtype.kind == "O":
-                # Something was added to the data array that caused numpy to want to
-                # treat it as a general-purpose Object array. This usually happens
-                # because:
-                # - the data contain NoDefaults
-                # - the data contain one or more Nones,
-                # - the data contain special types like tuples, dicts, etc
-                # - the data are composed of arrays that numpy would otherwise happily
-                # convert to a higher-order array, but the dimensions of the sub-arrays
-                # are inconsistent ("jagged")
-                # - there is some sort of honest-to-goodness weird object
-                # We want to support the first two cases with minimal intrusion, since
-                # these should be pretty easy to faithfully represent in the db. The
-                # jagged case should be supported as well, but may require a less
-                # faithful representation (e.g. flattened), but the last case isn't
-                # really worth supporting.
-
-                # Here is one proposal:
-                # - Check to see if the array is jagged. all(shape == shape[0]). If not,
-                # flatten, store the data offsets and array shapes, and None locations
-                # as attrs
-                # - If not jagged, all top-level ndarrays are the same shape, so it is
-                # easier to replace Nones with ndarrays filled with special values.
-                if parameters.NoDefault in data:
-                    data = None
-                else:
-                    data, specialAttrs = packSpecialData(data, paramDef.name)
-                    attrs.update(specialAttrs)
+            # - Check to see if the array is jagged. If so, flatten, store the
+            # data offsets and array shapes, and None locations as attrs.
+            # - If not jagged, all top-level ndarrays are the same shape, so it is
+            # easier to replace Nones with ndarrays filled with special values.
+            if isinstance(data, JaggedArray):
+                data, specialAttrs = packSpecialData(data, paramDef.name)
+                attrs.update(specialAttrs)
+
+            else:  # numpy.ndarray
+                # Convert Unicode to byte-string
+                if data.dtype.kind == "U":
+                    data = data.astype("S")
+
+                if data.dtype.kind == "O":
+                    # Something was added to the data array that caused numpy to want to
+                    # treat it as a general-purpose Object array. This usually happens
+                    # because:
+                    # - the data contain NoDefaults
+                    # - the data contain one or more Nones,
+                    # - the data contain special types like tuples, dicts, etc
+                    # - there is some sort of honest-to-goodness weird object
+                    # We want to support the first two cases with minimal intrusion, since
+                    # these should be pretty easy to faithfully represent in the db.
+                    # The last case isn't really worth supporting.
+
+                    if parameters.NoDefault in data:
+                        data = None
+                    else:
+                        data, specialAttrs = packSpecialData(data, paramDef.name)
+                        attrs.update(specialAttrs)
 
             if data is None:
                 continue
@@ -958,7 +969,9 @@ class Database3:
                         "should have been empty".format(paramDef.name, g)
                     )
 
-                dataset = g.create_dataset(paramDef.name, data=data, compression="gzip")
+                dataset = g.create_dataset(
+                    paramDef.name, data=data, compression="gzip", track_order=True
+                )
                 if any(attrs):
                     Database3._writeAttrs(dataset, h5group, attrs)
             except Exception:
@@ -982,7 +995,9 @@ class Database3:
         nDens = collectBlockNumberDensities(blocks)
 
         for nucName, numDens in nDens.items():
-            h5group.create_dataset(nucName, data=numDens, compression="gzip")
+            h5group.create_dataset(
+                nucName, data=numDens, compression="gzip", track_order=True
+            )
 
     @staticmethod
     def _readParams(h5group, compTypeName, comps, allowMissing=False):
@@ -1357,10 +1372,10 @@ class Database3:
                 # 3) not performing parameter renaming. This may become necessary
                 for paramName in params or h5GroupForType.keys():
                     if paramName == "location":
-                        # cast to a numpy array so that we can use list indices
-                        data = numpy.array(layout.location)[layoutIndicesForType][
-                            indexInData
-                        ]
+                        locs = []
+                        for id in indexInData:
+                            locs.append((layout.location[layoutIndicesForType[id]]))
+                        data = numpy.array(locs)
                     elif paramName in h5GroupForType:
                         dataSet = h5GroupForType[paramName]
                         try:
@@ -1396,7 +1411,9 @@ class Database3:
 
                     # iterating of numpy is not fast..
                     for c, val in zip(reorderedComps, data.tolist()):
-                        if isinstance(val, list):
+                        if paramName == "location":
+                            val = tuple(val)
+                        elif isinstance(val, list):
                             val = numpy.array(val)
 
                         histData[c][paramName][cycle, timeNode] = val
@@ -1444,7 +1461,7 @@ class Database3:
                 )
 
                 if "attrs" not in group:
-                    attrGroup = group.create_group("attrs")
+                    attrGroup = group.create_group("attrs", track_order=True)
                 else:
                     attrGroup = group["attrs"]
                 dataName = str(len(attrGroup)) + "_" + key
@@ -1490,7 +1507,7 @@ class Database3:
 
 
 def packSpecialData(
-    data: numpy.ndarray, paramName: str
+    arrayData: [numpy.ndarray, JaggedArray], paramName: str
 ) -> Tuple[Optional[numpy.ndarray], Dict[str, Any]]:
     """
     Reduce data that wouldn't otherwise play nicely with HDF5/numpy arrays to a format
@@ -1505,6 +1522,7 @@ def packSpecialData(
     composites that are storing the parameters. For instance, if we are dealing with a
     Block parameter, the first index in the numpy array of data is the block index; so
     if each block has a parameter that is a dictionary, ``data`` would be a ndarray,
+                    offset += arr.size
     where each element is a dictionary. This routine supports a number of different
     "strange" things:
 
@@ -1522,9 +1540,10 @@ def packSpecialData(
 
     Parameters
     ----------
-    data
-        An ndarray storing the data that we want to stuff into the database. These are
-        usually dtype=Object, which is how we usually end up here in the first place.
+    arrayData
+        An ndarray or JaggedArray object storing the data that we want to stuff into
+        the database. If the data is jagged, a special JaggedArray instance is passed
+        in, which contains a 1D array with offsets and shapes.
 
     paramName
         The parameter name that we are trying to store data for. This is mostly used for
@@ -1534,10 +1553,15 @@ def packSpecialData(
     --------
     unpackSpecialData
     """
-    # Check to make sure that we even need to do this. If the numpy data type is
-    # not "O", chances are we have nice, clean data.
-    if data.dtype != "O":
-        return data, {}
+    if isinstance(arrayData, JaggedArray):
+        data = arrayData.flattenedArray
+    else:
+        # Check to make sure that we even need to do this. If the numpy data type is
+        # not "O", chances are we have nice, clean data.
+        if arrayData.dtype != "O":
+            return arrayData, {}
+        else:
+            data = arrayData
 
     attrs: Dict[str, Any] = {"specialFormatting": True}
 
@@ -1560,7 +1584,7 @@ def packSpecialData(
     # A robust solution would need
     # to do this on a case-by-case basis, and re-do it any time we want to
     # write, since circumstances may change. Not only that, but we may need
-    # to do perform more that one of these operations to get to an array
+    # to perform more than one of these operations to get to an array
     # that we want to put in the database.
     if any(isinstance(d, dict) for d in data):
         # we're assuming that a dict is {str: float}. We store the union of
@@ -1585,6 +1609,13 @@ def packSpecialData(
 
         return data, attrs
 
+    if isinstance(arrayData, JaggedArray):
+        attrs["jagged"] = True
+        attrs["offsets"] = arrayData.offsets
+        attrs["shapes"] = arrayData.shapes
+        attrs["noneLocations"] = arrayData.nones
+        return data, attrs
+
     # conform non-numpy arrays to numpy
     for i, val in enumerate(data):
         if isinstance(val, (list, tuple)):
@@ -1595,45 +1626,6 @@ def packSpecialData(
         data = replaceNonesWithNonsense(data, paramName, nones)
         return data, attrs
 
-    # check if data is jagged
-    candidate = next((d for d in data if d is not None))
-    shape = candidate.shape
-    ndim = candidate.ndim
-    isJagged = (
-        not all(d.shape == shape for d in data if d is not None) or candidate.size == 0
-    )
-
-    if isJagged:
-        assert all(
-            val.ndim == ndim for val in data if val is not None
-        ), "Inconsistent dimensions in jagged array for: {}\nDimensions: {}".format(
-            paramName, [val.ndim for val in data if val is not None]
-        )
-        attrs["jagged"] = True
-
-        # offsets[i] is the index of the zero-th element of sub-array i
-        offsets = numpy.array(
-            [0]
-            + list(
-                itertools.accumulate(val.size if val is not None else 0 for val in data)
-            )[:-1]
-        )
-
-        # shapes[i] is the shape of the i-th sub-array. Nones are represented by all
-        # zeros
-        shapes = numpy.array(
-            list(val.shape if val is not None else ndim * (0,) for val in data)
-        )
-
-        data = numpy.delete(data, nones)
-
-        data = numpy.concatenate(data, axis=None)
-
-        attrs["offsets"] = offsets
-        attrs["shapes"] = shapes
-        attrs["noneLocations"] = nones
-        return data, attrs
-
     if any(isinstance(d, (tuple, list, numpy.ndarray)) for d in data):
         data = replaceNonesWithNonsense(data, paramName, nones)
         return data, attrs
@@ -1687,21 +1679,9 @@ def unpackSpecialData(data: numpy.ndarray, attrs, paramName: str) -> numpy.ndarr
     if attrs.get("jagged", False):
         offsets = attrs["offsets"]
         shapes = attrs["shapes"]
-        ndim = len(shapes[0])
-        emptyArray = numpy.ndarray(ndim * (0,), dtype=data.dtype)
-        unpackedJaggedData: List[Optional[numpy.ndarray]] = []
-        for offset, shape in zip(offsets, shapes):
-            if tuple(shape) == ndim * (0,):
-                # Start with an empty array. This may be replaced with a None later
-                unpackedJaggedData.append(emptyArray)
-            else:
-                unpackedJaggedData.append(
-                    numpy.ndarray(shape, dtype=data.dtype, buffer=data[offset:])
-                )
-        for i in attrs["noneLocations"]:
-            unpackedJaggedData[i] = None
-
-        return numpy.array(unpackedJaggedData, dtype=object)
+        nones = attrs["noneLocations"]
+        data = JaggedArray.fromH5(data, offsets, shapes, nones, data.dtype, paramName)
+        return data
     if attrs.get("dict", False):
         keys = numpy.char.decode(attrs["keys"])
         unpackedData = []
diff --git a/armi/bookkeeping/db/jaggedArray.py b/armi/bookkeeping/db/jaggedArray.py
new file mode 100644
index 00000000..3135f7da
--- /dev/null
+++ b/armi/bookkeeping/db/jaggedArray.py
@@ -0,0 +1,191 @@
+# Copyright 2024 TerraPower, LLC
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""
+Tooling to help flatten jagged (non rectangular) data into rectangular arrays.
+
+The goal here is to support jagged data for NumPy arrays to be written into the ARMI databases.
+"""
+
+from typing import List, Optional
+
+import numpy as np
+
+from armi import runLog
+
+
+class JaggedArray:
+    """
+    Take a list of numpy arrays or lists and flatten them into a single 1D array.
+
+    This implementation can preserve the structure of a multi-dimensional numpy array
+    by storing the dimensions in self.shapes and then re-populating a numpy array of
+    that shape from the flattened 1D array. However, it can only preserve one layer of
+    jaggedness in a list of lists (or other iterables). For example, a list of tuples
+    with varying lengths can be flattened and reconstituted exactly. But, if a list of
+    lists of tuples is passed in, the tuples in that final layer of nesting will all be
+    flattened to a single 1D numpy array after a round trip. No structure is retained
+    from nested lists of jagged lists or tuples.
+    """
+
+    def __init__(self, jaggedData, paramName):
+        """
+        JaggedArray constructor.
+
+        Parameters
+        ----------
+        jaggedData: list of np.ndarray
+            A list of numpy arrays (or lists or tuples) to be flattened into a single array
+        paramName: str
+            The name of the parameter represented by this data
+        """
+        offset = 0
+        flattenedArray = []
+        offsets = []
+        shapes = []
+        nones = []
+        for i, arr in enumerate(jaggedData):
+            if isinstance(arr, (np.ndarray, list, tuple)):
+                if len(arr) == 0:
+                    nones.append(i)
+                else:
+                    offsets.append(offset)
+                    try:
+                        numpyArray = np.array(arr)
+                        shapes.append(numpyArray.shape)
+                        offset += numpyArray.size
+                        flattenedArray.extend(numpyArray.flatten())
+                    except:  # noqa: E722
+                        # numpy might fail if it's jagged
+                        flattenedList = self.flatten(arr)
+                        shapes.append(
+                            len(flattenedList),
+                        )
+                        offset += len(flattenedList)
+                        flattenedArray.extend(flattenedList)
+            elif isinstance(arr, (int, float)):
+                offsets.append(offset)
+                shapes.append((1,))
+                offset += 1
+                flattenedArray.append(arr)
+            elif arr is None:
+                nones.append(i)
+
+        self.flattenedArray = np.array(flattenedArray)
+        self.offsets = np.array(offsets)
+        try:
+            self.shapes = np.array(shapes)
+        except ValueError as ee:
+            runLog.error(
+                "Error! It seems like ARMI may have tried to flatten a jagged array "
+                "where the elements have different numbers of dimensions. `shapes` "
+                "attribute of the JaggedArray for {} cannot be made into a numpy "
+                "array; it might be jagged.".format(paramName)
+            )
+            runLog.error(shapes)
+            raise ValueError(ee)
+        self.nones = np.array(nones)
+        self.dtype = self.flattenedArray.dtype
+        self.paramName = paramName
+
+    def __iter__(self):
+        """Iterate over the unpacked list."""
+        return iter(self.unpack())
+
+    def __contains__(self, other):
+        return other in self.flattenedArray
+
+    @staticmethod
+    def flatten(x):
+        """
+        Recursively flatten an iterable (list, tuple, or numpy.ndarray).
+
+        x : list, tuple, np.ndarray
+            An iterable. Can be a nested iterable in which the elements
+            themselves are also iterable.
+        """
+        if isinstance(x, (list, tuple, np.ndarray)):
+            if len(x) == 0:
+                return []
+            first, rest = x[0], x[1:]
+            return JaggedArray.flatten(first) + JaggedArray.flatten(rest)
+        else:
+            return [x]
+
+    @classmethod
+    def fromH5(cls, data, offsets, shapes, nones, dtype, paramName):
+        """
+        Create a JaggedArray instance from an HDF5 dataset.
+
+        The JaggedArray is stored in HDF5 as a flat 1D array with accompanying
+        attributes of "offsets" and "shapes" to define how to reconstitute the
+        original data.
+
+        Parameters
+        ----------
+        data: np.ndarray
+            A flattened 1D numpy array read in from an HDF5 file
+        offsets: np.ndarray
+            Offset indices for the zeroth element of each constituent array
+        shapes: np.ndarray
+            The shape of each constituent array
+        nones: np.ndarray
+            The location of Nones
+        dtype: np.dtype
+            The data type for the array
+        paramName: str
+            The name of the parameter represented by this data
+
+        Returns
+        -------
+        obj: JaggedArray An instance of JaggedArray populated with the input data
+        """
+        obj = cls([], paramName)
+        obj.flattenedArray = np.array(data)
+        obj.offsets = np.array(offsets)
+        obj.shapes = np.array(shapes)
+        obj.nones = np.array(nones)
+        obj.dtype = dtype
+        obj.paramName = paramName
+        return obj
+
+    def tolist(self):
+        """Alias for unpack() to make this class respond like a np.ndarray."""
+        return self.unpack()
+
+    def unpack(self):
+        """
+        Unpack a JaggedArray object into a list of arrays.
+
+        Returns
+        -------
+        unpackedJaggedData: list of np.ndarray
+            List of numpy arrays with varying dimensions (i.e., jagged arrays)
+        """
+        unpackedJaggedData: List[Optional[np.ndarray]] = []
+        numElements = len(self.offsets) + len(self.nones)
+        j = 0  # non-None element counter
+        for i in range(numElements):
+            if i in self.nones:
+                unpackedJaggedData.append(None)
+            else:
+                unpackedJaggedData.append(
+                    np.ndarray(
+                        self.shapes[j],
+                        dtype=self.dtype,
+                        buffer=self.flattenedArray[self.offsets[j] :],
+                    )
+                )
+                j += 1
+
+        return unpackedJaggedData
diff --git a/armi/physics/neutronics/__init__.py b/armi/physics/neutronics/__init__.py
index 463166c1..f995d488 100644
--- a/armi/physics/neutronics/__init__.py
+++ b/armi/physics/neutronics/__init__.py
@@ -36,8 +36,7 @@ from enum import IntEnum
 import numpy
 import tabulate
 
-from armi import plugins
-from armi import runLog
+from armi import plugins, runLog
 from armi.physics.neutronics.const import CONF_CROSS_SECTION
 
 
@@ -79,8 +78,8 @@ class NeutronicsPlugin(plugins.ArmiPlugin):
     @plugins.HOOKIMPL
     def defineSettings():
         """Define settings for the plugin."""
-        from armi.physics.neutronics import settings as neutronicsSettings
         from armi.physics.neutronics import crossSectionSettings
+        from armi.physics.neutronics import settings as neutronicsSettings
         from armi.physics.neutronics.fissionProductModel import (
             fissionProductModelSettings,
         )
@@ -99,10 +98,10 @@ class NeutronicsPlugin(plugins.ArmiPlugin):
     @plugins.HOOKIMPL
     def defineSettingsValidators(inspector):
         """Implementation of settings inspections for neutronics settings."""
-        from armi.physics.neutronics.settings import getNeutronicsSettingValidators
         from armi.physics.neutronics.fissionProductModel.fissionProductModelSettings import (
             getFissionProductModelSettingValidators,
         )
+        from armi.physics.neutronics.settings import getNeutronicsSettingValidators
 
         settingsValidators = getNeutronicsSettingValidators(inspector)
         settingsValidators.extend(getFissionProductModelSettingValidators(inspector))
@@ -133,7 +132,6 @@ from armi.physics.neutronics.const import (
     RESTARTFILES,
 )
 
-
 # ARC and CCCC cross section file format names
 COMPXS = "COMPXS"
 PMATRX = "PMATRX"
@@ -252,18 +250,20 @@ def applyEffectiveDelayedNeutronFractionToCore(core, cs):
         core.p.betaDecayConstants = numpy.array(decayConstants)
 
         reportTableData.append(("Total Delayed Neutron Fraction", core.p.beta))
-        reportTableData.append(
-            ("Group-wise Delayed Neutron Fractions", core.p.betaComponents)
-        )
-        reportTableData.append(
-            ("Group-wise Precursor Decay Constants", core.p.betaDecayConstants)
-        )
+        for i, betaComponent in enumerate(core.p.betaComponents):
+            reportTableData.append(
+                (f"Group {i} Delayed Neutron Fractions", betaComponent)
+            )
+        for i, decayConstant in enumerate(core.p.betaDecayConstants):
+            reportTableData.append(
+                ("Group {i} Precursor Decay Constants", decayConstant)
+            )
 
     # Report to the user the values were not applied.
     if not reportTableData and (beta is not None or decayConstants is not None):
         runLog.warning(
             f"Delayed neutron fraction(s) - {beta} and decay constants"
-            " - {decayConstants} have not been applied."
+            f" - {decayConstants} have not been applied."
         )
     else:
         runLog.extra(
diff --git a/armi/reactor/composites.py b/armi/reactor/composites.py
index 4ec1685e..d105cbd8 100644
--- a/armi/reactor/composites.py
+++ b/armi/reactor/composites.py
@@ -36,24 +36,19 @@ import collections
 import itertools
 import operator
 import timeit
-from typing import Dict, Optional, Type, Tuple, List, Union
+from typing import Dict, List, Optional, Tuple, Type, Union
 
 import numpy
 import six
 import tabulate
 
-from armi import context
-from armi import runLog
-from armi import utils
-from armi.nucDirectory import elements
-from armi.nucDirectory import nucDir, nuclideBases
+from armi import context, runLog, utils
+from armi.nucDirectory import elements, nucDir, nuclideBases
 from armi.physics.neutronics.fissionProductModel import fissionProductModel
-from armi.reactor import grids
-from armi.reactor import parameters
+from armi.reactor import grids, parameters
 from armi.reactor.flags import Flags, TypeSpec
 from armi.reactor.parameters import resolveCollections
-from armi.utils import densityTools
-from armi.utils import units
+from armi.utils import densityTools, units
 from armi.utils.densityTools import calculateNumberDensity
 
 
@@ -217,7 +212,7 @@ def _defineBaseParameters():
         the DB, it is possible to recover the flags from that.
         * Storing flags to the DB may be complicated, since it is easier to imagine a
         number of flags that is greater than the width of natively-supported integer
-        types, requiring some extra tricks to store the flages in an HDF5 file.
+        types, requiring some extra tricks to store the flags in an HDF5 file.
         * Allowing flags to be modified by plugins further complicates things, in that
         it is important to ensure that the meaning of all bits in the flag value are
         consistent between a database state and the current ARMI environment. This may
diff --git a/armi/utils/properties.py b/armi/utils/properties.py
index 2532b72e..6e7dd75d 100644
--- a/armi/utils/properties.py
+++ b/armi/utils/properties.py
@@ -28,8 +28,11 @@ def areEqual(val1, val2, relativeTolerance=0.0):
 def numpyHackForEqual(val1, val2):
     """Checks lots of types for equality like strings and dicts."""
     # when doing this with numpy arrays you get an array of booleans which causes the value error
-    notEqual = val1 != val2
+    if isinstance(val1, numpy.ndarray) and isinstance(val2, numpy.ndarray):
+        if val1.size != val2.size:
+            return False
 
+    notEqual = val1 != val2
     try:  # should work for everything but numpy arrays
         if isinstance(notEqual, numpy.ndarray) and notEqual.size == 0:
             return True
diff --git a/doc/release/0.3.rst b/doc/release/0.3.rst
index fa6d8054..86ffdfd8 100644
--- a/doc/release/0.3.rst
+++ b/doc/release/0.3.rst
@@ -14,6 +14,7 @@ New Features
 #. Reset ``Reactor`` data on worker processors after every interaction to free memory from state distribution.
    (`PR#1729 <https://github.com/terrapower/armi/pull/1729>`_ and `PR#1750 <https://github.com/terrapower/armi/pull/1750>`_)
 #. Density can be specified for components via ``custom isotopics`` in the blueprints. (`PR#1745 <https://github.com/terrapower/armi/pull/1745>`_)
+#. Implement a new `JaggedArray` class that handles HDF5 interface for jagged data. (`PR#1726 <https://github.com/terrapower/armi/pull/1726>`_)
 #. TBD
 
 API Changes
diff --git a/doc/user/outputs.rst b/doc/user/outputs.rst
index 8267fa93..24ec598b 100644
--- a/doc/user/outputs.rst
+++ b/doc/user/outputs.rst
@@ -248,5 +248,5 @@ still being able to faithfully reconstruct the original data. To accomplish this
 HDF5 dataset attributes to indicate when some manipulation is necessary. Writing
 such special data to the HDF5 file and reading it back again is accomplished with the
 :py:func:`armi.bookkeeping.db.database3.packSpecialData` and
-:py:func:`armi.bookkeeping.db.database3.packSpecialData`. Refer to their implementations
+:py:func:`armi.bookkeeping.db.database3.unpackSpecialData`. Refer to their implementations
 and documentation for more details.
diff --git a/pyproject.toml b/pyproject.toml
index 6472e7ed..3aa0251f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,7 +34,7 @@ dependencies = [
     "h5py>=3.0,<=3.9", # Needed because our database files are H5 format
     "htmltree>=0.7.6", # Our reports have HTML output
     "matplotlib>=3.5.3,<3.8.0", # Important plotting library
-    "numpy>=1.21,<=1.23.5", # Important math library
+    "numpy>=1.21", # Important math library
     "ordered-set>=3.1.1", # A useful data structure
     "pluggy>=1.2.0", # Central tool behind the ARMI Plugin system
     "pyDOE>=0.3.8", # We import a Latin-hypercube algorithm to explore a phase space


Test patch:
diff --git a/armi/bookkeeping/db/tests/test_database3.py b/armi/bookkeeping/db/tests/test_database3.py
index 3efe91c8..18b8921b 100644
--- a/armi/bookkeeping/db/tests/test_database3.py
+++ b/armi/bookkeeping/db/tests/test_database3.py
@@ -21,6 +21,7 @@ import numpy
 
 from armi.bookkeeping.db import _getH5File
 from armi.bookkeeping.db import database3
+from armi.bookkeeping.db.jaggedArray import JaggedArray
 from armi.bookkeeping.db.databaseInterface import DatabaseInterface
 from armi.reactor import parameters
 from armi.reactor.tests.test_reactors import loadTestReactor, reduceTestReactorRings
@@ -324,13 +325,11 @@ class TestDatabase3Smaller(unittest.TestCase):
         data2fNones = numpy.array(
             [None, [[1.0, 2.0, 6.0], [2.0, 3.0, 4.0]]], dtype=object
         )
-        dataJag = numpy.array(
-            [[[1, 2], [3, 4]], [[1, 2, 3], [4, 5, 6], [7, 8, 9]]], dtype=object
-        )
-        dataJagNones = numpy.array(
-            [[[1, 2], [3, 4]], [[1], [1]], [[1, 2, 3], [4, 5, 6], [7, 8, 9]]],
-            dtype=object,
-        )
+        twoByTwo = numpy.array([[1, 2], [3, 4]])
+        twoByOne = numpy.array([[1], [None]])
+        threeByThree = numpy.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+        dataJag = JaggedArray([twoByTwo, threeByThree], "testParam")
+        dataJagNones = JaggedArray([twoByTwo, twoByOne, threeByThree], "testParam")
         dataDict = numpy.array(
             [{"bar": 2, "baz": 3}, {"foo": 4, "baz": 6}, {"foo": 7, "bar": 8}]
         )
diff --git a/armi/bookkeeping/db/tests/test_jaggedArray.py b/armi/bookkeeping/db/tests/test_jaggedArray.py
new file mode 100644
index 00000000..89d6c53f
--- /dev/null
+++ b/armi/bookkeeping/db/tests/test_jaggedArray.py
@@ -0,0 +1,124 @@
+# Copyright 2019 TerraPower, LLC
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+"""Tests for the JaggedArray class."""
+import unittest
+
+import h5py
+import numpy
+
+from armi.bookkeeping.db.jaggedArray import JaggedArray
+from armi.utils.directoryChangers import TemporaryDirectoryChanger
+
+
+class TestJaggedArray(unittest.TestCase):
+    """Tests for the JaggedArray class."""
+
+    def setUp(self):
+        self.td = TemporaryDirectoryChanger()
+        self.td.__enter__()
+
+    def tearDown(self):
+        self.td.__exit__(None, None, None)
+
+    def test_roundTrip(self):
+        """Basic test that we handle Nones correctly in database read/writes."""
+        dataSet = [1, 2.0, None, [], [3, 4], (5, 6, 7), numpy.array([8, 9, 10, 11])]
+        self._compareRoundTrip(dataSet, "test-numbers")
+
+    def test_roundTripBool(self):
+        """Basic test that we handle Nones correctly in database read/writes."""
+        dataSet = [True, True, [False, True, False]]
+        self._compareRoundTrip(dataSet, "test-bool")
+
+    def test_flatten(self):
+        """Test the recursive flattening static method."""
+        testdata = [(1, 2), [3, 4, 5], [], None, 6, numpy.array([7, 8, 9])]
+        flatArray = JaggedArray.flatten(testdata)
+        self.assertEqual(flatArray, [1, 2, 3, 4, 5, None, 6, 7, 8, 9])
+
+    def _compareRoundTrip(self, data, paramName):
+        """Make sure that data is unchanged by packing/unpacking."""
+        jaggedArray = JaggedArray(data, paramName)
+
+        # write to HDF5
+        h5file = "test_jaggedArray.h5"
+        with h5py.File(h5file, "w") as hf:
+            dset = hf.create_dataset(
+                data=jaggedArray.flattenedArray,
+                name=jaggedArray.paramName,
+            )
+            dset.attrs["jagged"] = True
+            dset.attrs["offsets"] = jaggedArray.offsets
+            dset.attrs["shapes"] = jaggedArray.shapes
+            dset.attrs["noneLocations"] = jaggedArray.nones
+
+        with h5py.File(h5file, "r") as hf:
+            dataset = hf[paramName]
+            values = dataset[()]
+            offsets = dataset.attrs["offsets"]
+            shapes = dataset.attrs["shapes"]
+            nones = dataset.attrs["noneLocations"]
+
+        roundTrip = JaggedArray.fromH5(
+            values,
+            offsets,
+            shapes,
+            nones,
+            dtype=jaggedArray.flattenedArray.dtype,
+            paramName=paramName,
+        )
+        self._compareArrays(data, roundTrip)
+
+    def _compareArrays(self, ref, src):
+        """
+        Compare two numpy arrays.
+
+        Comparing numpy arrays that may have unsavory data (NaNs, Nones, jagged
+        data, etc.) is really difficult. For now, convert to a list and compare
+        element-by-element.
+
+        Several types of data do not survive a round trip. The if-elif branch
+        here converts the initial data into the format expected to be produced
+        by the round trip. The conversions are:
+
+        - For scalar values (int, float, etc.), the data becomes a numpy
+          array with a dimension of 1 after the round trip.
+        - Tuples and lists become numpy arrays
+        - Empty lists become `None`
+
+        """
+        # self.assertEqual(type(src), JaggedArray)
+        if isinstance(ref, numpy.ndarray):
+            ref = ref.tolist()
+            src = src.tolist()
+
+        for v1, v2 in zip(ref, src):
+            # Entries may be None
+            if isinstance(v1, numpy.ndarray):
+                v1 = v1.tolist()
+            elif isinstance(v1, tuple):
+                v1 = list(v1)
+            elif isinstance(v1, int):
+                v1 = numpy.array([v1])
+            elif isinstance(v1, float):
+                v1 = numpy.array([v1], dtype=numpy.float64)
+            elif v1 is None:
+                pass
+            elif len(v1) == 0:
+                v1 = None
+
+            if isinstance(v2, numpy.ndarray):
+                v2 = v2.tolist()
+
+            self.assertEqual(v1, v2)
diff --git a/armi/nuclearDataIO/tests/test_xsLibraries.py b/armi/nuclearDataIO/tests/test_xsLibraries.py
index 35d7f606..3fe2bc81 100644
--- a/armi/nuclearDataIO/tests/test_xsLibraries.py
+++ b/armi/nuclearDataIO/tests/test_xsLibraries.py
@@ -18,13 +18,12 @@ import os
 import traceback
 import unittest
 
+import numpy
 from six.moves import cPickle
 
 from armi.nucDirectory import nuclideBases
 from armi.nuclearDataIO import xsLibraries
-from armi.nuclearDataIO.cccc import gamiso
-from armi.nuclearDataIO.cccc import isotxs
-from armi.nuclearDataIO.cccc import pmatrx
+from armi.nuclearDataIO.cccc import gamiso, isotxs, pmatrx
 from armi.tests import mockRunLogs
 from armi.utils import properties
 from armi.utils.directoryChangers import TemporaryDirectoryChanger
@@ -356,8 +355,8 @@ class TestXSlibraryMerging(TempFileMixin):
 
     def test_cannotMergeXSLibxWithDifferentGroupStructure(self):
         dummyXsLib = xsLibraries.IsotxsLibrary()
-        dummyXsLib.neutronEnergyUpperBounds = [1, 2, 3]
-        dummyXsLib.gammaEnergyUpperBounds = [1, 2, 3]
+        dummyXsLib.neutronEnergyUpperBounds = numpy.array([1, 2, 3])
+        dummyXsLib.gammaEnergyUpperBounds = numpy.array([1, 2, 3])
         with self.assertRaises(properties.ImmutablePropertyError):
             dummyXsLib.merge(self.libCombined)
 
@@ -439,7 +438,7 @@ class Pmatrx_merge_Tests(TestXSlibraryMerging):
 
     def test_cannotMergeXSLibsWithDifferentGammaGroupStructures(self):
         dummyXsLib = xsLibraries.IsotxsLibrary()
-        dummyXsLib.gammaEnergyUpperBounds = [1, 2, 3]
+        dummyXsLib.gammaEnergyUpperBounds = numpy.array([1, 2, 3])
         with self.assertRaises(properties.ImmutablePropertyError):
             dummyXsLib.merge(self.libCombined)