.. ======================================================================
..  © (or copyright) 2019-2026. Triad National Security, LLC. All rights
..  reserved.  This program was produced under U.S. Government contract
..  89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
..  operated by Triad National Security, LLC for the U.S.  Department of
..  Energy/National Nuclear Security Administration. All rights in the
..  program are reserved by Triad National Security, LLC, and the
..  U.S. Department of Energy/National Nuclear Security
..  Administration. The Government is granted for itself and others acting
..  on its behalf a nonexclusive, paid-up, irrevocable worldwide license
..  in this material to reproduce, prepare derivative works, distribute
..  copies to the public, perform publicly and display publicly, and to
..  permit others to do so.
.. ======================================================================

.. Generative AI was used to assist with modifications to this document.

.. _databox:

The DataBox
===========

The fundamental data type in ``spiner`` is the ``DataBox``. A
``DataBox`` packages a multi-dimensional (up to six dimensions) array
with routines for interpolating on the array and for saving the data
to and loading the data from file.

To use databox, simply include the relevant header:

.. code-block:: cpp

  #include <databox.hpp>

``DatBox`` is templated on underyling data type, which defaults to the
``Real`` type provided by ``ports-of-call``. (This is usually a
``double``.)

.. note::
  The default type can be set to type ``float`` if the preprocessor
  macro ``SINGLE_PRECISION_ENABLED`` is defined.

Any arithmetic type is supported, although the code has only been
tested carefully with floating point numbers. To set ``DataBox`` to a
single type, you may wish to declare a type alias such as:

.. code-block:: cpp

   using DataBox = Spiner::DataBox<double>

Spiner is also templated on how the interpolation gridding works. This
template parameter is called ``Grid_t``. The available options at this time are:

* ``Spiner::RegularGrid1D<T>``
* ``Spiner::PiecewiseGrid1D<T>``

where here ``T`` is the arithmetic type as discussed above. The
default type is ``RegularGrid1D``. You can further alias ``DataBox``
as, for example:

.. code-block:: cpp

   using DataBox = Spiner::DataBox<double, Spiner::RegularGrid1D<double>>;

More detail on the interpolation gridding is available below and in
the interpolation section.

The third template parameter, ``GradStencil``, selects the implementation of
``interpToRealWithGrads`` and defaults to
``Spiner::InterpolationStencils::Linear``:

.. code-block:: cpp

   using DataBox = Spiner::DataBox<double, Spiner::RegularGrid1D<double>,
                                  MyGradStencil>;

The built-in ``Spiner::InterpolationStencils::Rational`` policy enables
:ref:`rational interpolation <rational-interpolation>` with gradients in
1D and 2D.

The policy supplies portable static overloads
``interpToRealWithGrads(x1, grids, data)`` and
``interpToRealWithGrads(x2, x1, grids, data)``, returning
``Spiner::InterpolationHelpers::Result1D<T>`` and ``Result2D<T>`` respectively.
``grids`` and ``data`` are passed by const reference; grid index zero is the
fastest axis, and data are accessed as ``data(i1)`` or ``data(i2, i1)``.
Policies must support the calling execution space and must not throw, since
the DataBox entry points are ``noexcept``. DataBox retains the rank and grid
validity assertions before dispatching to the policy. The policy type is
available as ``DataBox::GradStencilType`` and is preserved by slices, copies,
and device transfers. It does not change value-only ``interpToReal`` calls.
The internal ``Concept`` template parameter follows ``GradStencil`` as the
fourth parameter; code explicitly supplying the former third parameter must
move it to fourth position.

.. note::
   In C++17 and later, you can also get the default type specialization
   by simply omitting the template arguments.

.. note::
  In the function signatures below, GPU/performance portability
  decorators have been excluded for brevity. However they are present
  in the actual code.

.. note::
   In the function signatures below, we will often refer to the type
   ``Real`` and the type ``T``. These are both references to the
   underlying templated arithmetic type.
   
Creating a ``DataBox``
^^^^^^^^^^^^^^^^^^^^^^

You can create a ``DataBox`` of a given shape via the constructor:

.. code-block:: cpp

  int nx1 = 2;
  int nx2 = 3;
  int nx3 = 4;
  Spiner::DataBox<double> db(nx3, nx2, nx1);

The constructor takes any number of shape values (e.g., ``nx*``) up to
six (or ``Spiner::MAXRANK``) values. Zero shape values initializes an
empty, size-zero array.

.. note::
  ``DataBox`` is row-major ordered. By convention, ``x3`` is the
  slowest moving index and ``x1`` is the fastest.

If GPU support is enabled, a ``DataBox`` can be allocated on either
host or device, depending on the ``AllocationTarget``. For example, to
explicitly allocate one array on the host and one on the device, you
might call:

.. code-block:: cpp

  // Allocates on the host (CPU)
  Spiner::DataBox<double> db_host(Spiner::AllocationTarget::Host, nx2, nx1);
  // Allocates on the device (GPU)
  Spiner::DataBox<double> db_dev(Spiner::AllocationTarget::Device, nx2, nx1);

.. note::
  If GPU support is not enabled, these both allocate on host.

You can also wrap a ``DataBox`` around a pointer you allocated
yourself. For example:

.. code-block:: cpp

  std::vector<double> mydata(nx1*nx2);
  Spiner::DataBox<double> db(mydata.data(), nx2, nx1);

You can also resize a ``DataBox``, which you can use to modify a
``DataBox`` in-place. For example:

.. code-block:: cpp

  Spiner::DataBox<double> db; // empty
  // clears old memory, resizes the underlying array,
  // and resets strides
  db.resize(nx3, nx2, nx1);
 
Just like the constructor, ``resize`` takes an optional (first)
argument for the ``AllocationTarget``.

.. warning::
  ``DataBox::resize`` is destructive. The underlying data is not preserved.

If you want to change the stride without changing the underlying data,
you can use ``reshape``, which modifies the dimensions of the
array, without modifying the underlying memory. For example:

.. code-block:: cpp

  // allocate a 1D databox
  Spiner::DataBox<double> db(nx3*nx2*nx1);
  // interpret it as a 3D object
  db.reshape(nx3, nx2, nx1);

.. warning::

  Make sure not to change the underlying size of the array
  when using ``reshape``. This is checked with an ``assert``
  statement, so you will get errors when compiling without
  the ``NDEBUG`` preprocessor macro.

The method

.. cpp:function:: void DataBox::reset();

sets the ``DataBox`` to be empty with zero rank.

Copying a ``DataBox`` to device
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

If GPU support is enabled, you can deep-copy a ``DataBox``, its tabulated
values, and (by default) its grids from host to device with the function

.. cpp:function:: DataBox getOnDeviceDataBox(const DataBox &db_host, bool include_grids = true);

which returns a new databox with the data in ``db_host`` copied to
GPU. An object-oriented method

.. cpp:function:: DataBox DataBox::getOnDevice(bool include_grids = true) const;

exists as well, which returns a new object with the underlying data
copied to GPU.

Set ``include_grids`` to ``false`` to skip the deep copy of the grid
objects. This permits manual management of grid memory when reducing the
device-memory footprint is important. The result holds non-owning shallow
copies of the source's grids (see ``shallowCopy`` in the grid
documentation), so finalizing them does not free the source's grid memory.
Its grids must be made device-accessible before they are used for
device-side interpolation.

.. note::
  If GPU support is not enabled, ``getOnDevice`` and friends are
  no-ops.

Semantics and Memory Management
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

``DataBox`` has reference semantics---meaning that copying a
``DataBox`` does not copy the underlying data. In other words,

.. code-block:: cpp

  Spiner::DataBox<double> db1(size);
  Spiner::DataBox<double> db2 = db1;

shallow-copies ``db1`` into ``db2``. Especially for `Kokkos`_ like
workflows, this is very useful.

.. _`Kokkos`: https://github.com/kokkos/kokkos

.. warning::
  ``DataBox`` is neither reference-counted nor garbage-collected.
  If you create a ``DataBox`` you must clear the memory allocated
  just like you would for a pointer.

Two functions are provided for freeing memory in ``DataBox``:

.. cpp:function:: void free(DataBox& db);

and

.. cpp:function:: DataBox::finalize();

both will do the same thing and free the values in a ``DataBox`` in a
context-dependent way. I.e., no matter what the ``AllocationTarget``
was, the appropriate memory will be freed.

Some grids, such as ``NonUniformGrid1D``, own dynamic memory of their
own. That memory is managed separately from the values and is **not**
freed by ``free`` or ``finalize``. Free it with

.. cpp:function:: void DataBox::finalizeGrids();

which finalizes every grid held by the ``DataBox``. Grids that do not
own memory, including all ``RegularGrid1D`` and ``PiecewiseGrid1D``
grids and non-owning shallow copies, are left alone, so it is always
safe to call. It is independent of ``finalize`` and may be called on a
``DataBox`` that does not own its values. For example:

.. code-block:: cpp

  Spiner::DataBox<double, Spiner::NonUniformGrid1D<double>> db(N);
  db.setRange(0, points); // db's grid owns a copy of points
  // ...
  db.finalizeGrids(); // frees the grid coordinates
  free(db);           // frees the values

Grid ownership follows one rule. An operation that makes a new handle
to the *same* ``DataBox`` copies its grids along with their ownership.
An operation that builds a *different* ``DataBox`` from another's
metadata either gives it non-owning shallow copies of the grids or, if
it is a deep copy, new grids that it owns. It never makes two
``DataBox`` objects owners of the same grid memory. Specifically:

* ``setRange`` copies the grid it is given, including its ownership,
  so a ``DataBox`` given a memory-owning grid is responsible for it.
* Copy construction and assignment of a ``DataBox`` produce handles to
  the same ``DataBox``. They report the same ownership as the source,
  for both values and grids. Free through only one of them.
* Slices, ``copyShape``, and ``getOnDevice(false)`` build a different
  ``DataBox`` and give it non-owning shallow copies of the grids.
* ``getOnDevice()``, ``copyMetadata``, and ``deepCopy`` build a
  different ``DataBox`` with new grids that it owns.

.. warning::
  Do not free a ``DataBox`` if its memory is managed externally, e.g.,
  via a ``std::vector``. ``DataBox`` checks for this use-case
  via an ``assert`` statement.

You can check whether a given ``DataBox`` is empty, unmanaged, or
allocated on host or device with the

.. cpp:function:: DataBox::dataStatus() const;

method. It returns an ``enum class``, ``Spiner::DataStatus``, which
can take on the values ``Empty``, ``Unmanaged``, ``AllocatedHost``, or
``AllocatedDevice``. You can also check whether or not ``free`` should
be called with the method

.. cpp:function:: bool DataBox::ownsAllocatedMemory();

which returns ``true`` if a given databox is managing memory and
``false`` otherwise. The method

.. cpp:function:: bool DataBox::isReference();

returns ``false`` if the databox is managing memory and ``true``
otherwise.

To make an independent copy of a host-accessible ``DataBox``, use

.. cpp:function:: void DataBox::deepCopy(const DataBox &src);

which allocates new memory and copies the values, index types, and
grids of ``src``. The grids are deep copies owned by the new
``DataBox``, so free it with both ``finalizeGrids`` and ``finalize``.
The grids of the ``DataBox`` being copied into must not own memory.
Use ``getOnDevice`` to copy a ``DataBox`` to device.

Using ``DataBox`` with smart pointers
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Smart pointers can be used to manage a ``DataBox`` and automatically
call ``free`` for you, so long as you use them with a custom
deleter. Spiner provides the following deleter for use in this
scenario:

.. code-block:: cpp

  template <bool FinalizeGrids = false>
  struct DBDeleter {
    template <typename T>
    void operator()(T *ptr) {
      if constexpr (FinalizeGrids) {
        ptr->finalizeGrids();
      }
      ptr->finalize();
      delete ptr;
    }
  };

``DBDeleter<>`` frees only the values, like ``free``.
``DBDeleter<true>`` also calls ``finalizeGrids``, which you need when
the ``DataBox`` owns memory-owning grids.

It can be used, for example, with a ``std::unique_ptr`` via:

.. code-block:: cpp

  // needed for smart pointers
  #include <memory>

  // Creates a unique pointer pointing to a DataBox
  // with memory allocated on device
  std::unique_ptr<DataBox, Spiner::DBDeleter<>> pdb(
    new DataBox(Spiner::AllocationTarget::Device, N));
  
  // Before using the databox in, e.g., a GPU or Kokkos kernel, get a
  // shallow copy:
  auto db = *pdb;
  // some kokkos code...
  
  // when you leave scope, the data box will be freed.

Serialization and de-serialization
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Shared memory models, such as `MPI Windows`_, require allocation of
memory through an external API call (e.g.,
``MPI_Win_allocate_shared``), which tabulated data must be written
to. ``Spiner`` supports this model through **serialization** and
**de-serialization**. The relevant methods are as follows. The
function

.. cpp:function:: std::size_t DataBox::serializedSizeInBytes() const;

reports how much memory a ``DataBox`` object requires to be externally
allocated. The function

.. cpp:function:: std::size_t DataBox::dynamicMemorySizeInBytes() const;

reports the size of the value buffer and grid-owned dynamic memory,
excluding the inline ``DataBox`` object representation. The function

.. cpp:function:: std::size_t DataBox::dumpDynamicMemory(std::byte *dst) const;

writes only that dynamic memory. It is used internally by ``serialize``
after the inline object bytes have been copied.

.. cpp:function:: std::size_t serialize(std::byte *dst) const;

takes a ``std::byte*`` pointer, assumed to contain enough space for a
``DataBox``, and stores all information needed for the ``DataBox`` to
reconstruct itself. The return value is the amount of memory in bytes
used in the array by the serialized ``DataBox`` object. This method is
non-destructive; the original ``DataBox`` is unchanged. The function

.. cpp:function:: std::size_t DataBox::setPointer(T *src);

with the overload

.. cpp:function:: std::size_t DataBox::setPointer(std::byte *src);

sets the underlying tabulated data from the src pointer, which is
assumed to be the right size and shape, and relocates the trailing
grid-owned dynamic memory. This is useful for the deSerialize
function (described below) and for building your own
serialization/de-serialization routines in composite objects. The function

.. cpp:function:: std::size_t DataBox::deSerialize(std::byte *src);

initializes a ``DataBox`` to match the serialized ``DataBox``
contained in the ``src`` pointer.

.. note::

  Note that the de-serialized ``DataBox`` has **unmanaged** memory, as
  it is assumed that the ``src`` pointer manages its memory for
  it. Therefore, one **cannot** ``free`` the ``src`` pointer until
  everything you want to do with the de-serialized ``DataBox`` is
  over.

The serialized memory is laid out as the inline ``DataBox`` metadata,
followed by the tabulated values and grid-owned dynamic memory. The
inline grid objects are already part of the ``DataBox`` metadata, so
their static bytes are not emitted a second time. Dynamic memory is
dumped for every dimension, including indexed and named dimensions, so
resource management does not depend on mutable interpolation metadata.

The source buffer must remain alive and suitably aligned for the
entire lifetime of every ``DataBox`` de-serialized from it.

.. warning::

  ``DataBox`` previously exposed the serialize and deserialize
  routines using the C-style ``char*`` pointer. ``spiner`` now prefers
  the more ``std::byte*`` pointer. The former is still supported via
  overloads but should be considered deprecated.

Putting this all together, an application of
serialization/de-serialization probably looks like this:

.. code-block:: cpp

  // load a databox from, e.g., file
  Spiner::DataBox<double> db;
  db.loadHDF(filename);
  
  // get size of databox
  std::size_t allocate_size = db.serializedSizeInBytes();
  
  // Allocate the memory for the new databox.
  // In practice this would be an API call for, e.g., shared memory
  std::byte *memory = (std::byte*)malloc(allocate_size);
  
  // serialize the old databox
  std::size_t write_size = db.serialize(memory);
  
  // make a new databox and de-serialize it
  Spiner::DataBox<double> db2;
  std::size_t read_size = db2.deSerialize(memory);

  // read_size, write_size, and allocate_size should all be the same.
  assert((read_size == write_size) && (write_size == allocate_size));

.. warning::

  ``DataBox`` serialization is intended for transient communication
  between processes using the same Spiner version and a compatible
  architecture. The binary representation is not guaranteed to be
  compatible across Spiner versions, architectures, compilers, or
  build configurations. This is due to, for example, differences in
  endianness, alignment, and padding.

  Do not use serialization as a persistent-storage or web-interchange
  format. Use HDF5 for persistent, portable storage. Size calculation,
  serialization, and de-serialization must all use a compatible Spiner
  build.

.. _`MPI Windows`: https://www.mpi-forum.org/docs/mpi-4.1/mpi41-report/node311.htm

Accessing Elements of a ``DataBox``
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Elements of a ``DataBox`` can be accessed and set via the ``()``
operator. For example:

.. code-block:: cpp

  Spiner::DataBox<double> db(nx3, nx2, nx1);
  db(2,1,0) = 5.0;

The ``()`` operator accepts between one and six indexes. If you pass
in more indexes than the rank of the array, the excess indices are
ignored. If you pass in fewer, the unset indices are assumed to be
zero. The exception is the one-dimensional operator. You can always
stride through the "flattened" array by using the one-dimensional
accessor. For example:

.. code-block:: cpp

  for (int i = 0; i < nx3*nx2*nx1; ++i) {
    db(i) = static_cast<double>(i);
  }

fills the three-dimensional array above with the flat index of each
element.

Slicing
^^^^^^^^

A new ``DataBox`` containing a shallow slice of another ``DataBox``
can be constructed with the ``slice`` method:

.. cpp:function:: DataBox DataBox::slice(const int dim, const int indx, const int nvar) const;

this is fairly limited functionality. It returns a new ``DataBox``
containing only elements from ``indx`` to ``indx + nvar - 1`` in the
``dim`` direction. All other directions are unchanged. The slowest
moving dimension can be sliced to a single index with

.. cpp:function:: DataBox DataBox::slice(const int indx) const;
   
and the slowst-moving two dimensions can be sliced to a single pair of
indicies with

.. cpp:function:: DataBox DataBox::slice(const int i2, int i1) const;

Index Types and Interpolation Ranges
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Often-times an array mixes "continuous" and "discrete" variables. In
other words, some indices of an array are discretizations of a
continuous quantity, and we want to interpolate in those directions,
but other indices are discrete---they may index a particle species,
for example. A common example is in neutrino transport, where an array
of emissivities may depend on fluid density, fluid temperature,
electron fraction, neutrino energy, and neutrino species. The species
can only take three discrete values, but the density, temperature, and
electron fraction are all continuous.

``Spiner`` accounts for this by assigning each dimension in the array
a "type," represented as an ``enum class``, ``IndexType``. Currently
the type can be either ``Interpolated`` or ``Indexed``. When a new
``DataBox`` is created, all dimensions are set to
``IndexType::Indexed``. A dimension can be set to ``Interpolated`` via
the ``setRange`` method.

.. cpp:function:: void DataBox::setRange(int i, Grid_t g);

where here ``i`` is the dimension and ``g`` is the gridding object for
this index. In the default setup, where grids are uniformly spaced
(i.e., you use a ``RegularGrid1D``), this is:

.. cpp:function:: void DataBox::setRange(int i, T min, T max, int N);
   
where here ``i`` is the dimension, ``min`` is the minimum value of the
independent variable, ``max`` is the maximum value of the indpendent
variable, and ``N`` is the number of points in the ``i``
dimension. (Here ``T`` is the underlying templated data type.)

.. warning::
  In these routines, the dimension is indexed from zero. Moreover, the
  indexing possibly the opposite direction you're expecting. The zeroth
  dimension is the right most index when ``interpToReal`` is called.

.. note::
  There is a set of lower-level objects, ``RegularGrid1D``, and
  ``PiecewiseGrid1D``, which represent these interpolation ranges
  internally. There is a getter method ``range`` that works
  with the underlying ``Grid_t`` class directly. For
  more details, see the relevant documentation.

It's often desirable to have multiple databoxes with the exact same
shape and interpolation structure (i.e., independent variable
ranges). In this case, the method

.. cpp:function:: void DataBox::copyMetadata(const DataBox &src);

can assist. This method resets and re-allocates the data in a
``DataBox`` to the exact same size and shape as ``src``. More
importantly, it also copies the relevant ``IndexType`` and independent
variable range for each dimension. The grids are deep copies owned by
the ``DataBox``, so free them with ``finalizeGrids``. The grids of the
``DataBox`` being reset must not own memory.

One can also manually set the ``IndexType`` in a given dimension with

.. cpp:function:: void DataBox::setIndexType(int i, IndexType t);
   
and retrieve the ``IndexType`` with

.. cpp:function:: IndexType &DataBox::indexType(const int i);

to see if a dimension is interpolatable.

Interpolation to a real number
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The family of ``DataBox::interpToReal`` methods interpolate the
"entire" ``DataBox`` to a real number. Up to four-dimensional
interpolation is supported:

.. cpp:function:: T DataBox::interpToReal(const T x) const;

.. cpp:function:: T DataBox::interpToReal(const T x2, const T x1) const;

.. cpp:function:: T DataBox::interpToReal(const T x3, const T x2, const T x1) const;

.. cpp:function:: T DataBox::interpToReal(const T x4, const T x3, const T x2, const T x1) const;

where ``x1`` is the fastest moving direction, ``x2`` is less fast, and
so on. These interpolation routines are hand-tuned for performance.

.. warning::
  Do not call ``interpToReal`` with a ``DataBox`` that is the wrong shape
  or try to interpolate on indices that are not interpolatable.
  This is checked with an ``assert`` statement.

.. _interpolation-with-gradients:

Interpolation with gradients
^^^^^^^^^^^^^^^^^^^^^^^^^^^^

``interpToRealWithGrads`` returns the interpolated value and its derivatives
with respect to the continuous grid coordinates. It supports fully
interpolated rank-1 and rank-2 databoxes:

.. cpp:function:: DataBox::InterpResult1D DataBox::interpToRealWithGrads(const T x1) const noexcept;

.. cpp:function:: DataBox::InterpResult2D DataBox::interpToRealWithGrads(const T x2, const T x1) const noexcept;

The nested result aliases name plain aggregates with members of type ``T``:

.. code-block:: cpp

   struct InterpResult1D { T value, d_dx1; };
   struct InterpResult2D { T value, d_dx2, d_dx1; };

Use named members or structured bindings. As with all databox types, arguments
and returned derivatives follow slowest-to-fastest order: ``x2`` then
``x1``. For example, this host-side table represents
``f(y, x) = 1 + 2*x - 3*y + 4*x*y``:

.. code-block:: cpp

   Spiner::DataBox<double> db(3, 4); // (ny, nx)
   db.setRange(0, 0.0, 3.0, 4);    // x, fastest dimension
   db.setRange(1, 0.0, 4.0, 3);    // y, slowest dimension
   for (int j = 0; j < 3; ++j) {
     const double y = 2.0 * j;
     for (int i = 0; i < 4; ++i) {
       const double x = i;
       db(j, i) = 1.0 + 2.0*x - 3.0*y + 4.0*x*y;
     }
   }
   auto [value, df_dy, df_dx] = db.interpToRealWithGrads(1.0, 0.5);
   // value = 1, df_dy = -1, df_dx = 6
   db.finalize();

For a fully interpolated 1D databox, use:

.. code-block:: cpp

   auto [value, df_dx] = db1.interpToRealWithGrads(x);
   // Equivalently:
   auto result = db1.interpToRealWithGrads(x);
   // result.value and result.d_dx1

The derivatives are analytic derivatives of the selected policy's
interpolant. With the default ``Linear`` policy, this uses the same cell
selection and weights as ``interpToReal``. They are not interpolated samples of a separately
computed gradient field.  At knots or cell boundaries, derivatives can
be discontinuous; the result uses the derivative from the cell
selected by the grid lookup. With ``Linear``, out-of-range coordinates use the boundary
cell's linear extrapolation and its derivatives.

Derivatives are taken with respect to the supplied grid coordinates. If a
coordinate represents a transformed variable, such as a logarithm, apply
the chain rule to obtain derivatives with respect to the original variable.

All four built-in grid types support these methods. Custom grid types
must provide :ref:`weightsWithGrad <grid-weight-derivatives>` to
support this mode for linear interpolants. The databox must have the
matching rank and all dimensions must be interpolatable with
well-formed grids, as checked by assertions. There are no mixed
indexed/interpolated or higher-dimensional overloads of
``interpToRealWithGrads``.

.. _rational-interpolation:

Rational function interpolation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

``Spiner::InterpolationStencils::Rational`` implements the rational function
interpolation described by
`Kerley (1977, LA-6903-MS) <https://www.osti.gov/servlets/purl/7301808>`_,
with analytic coordinate derivatives. Select it as the third ``DataBox``
template parameter:

.. code-block:: cpp

   using RationalBox = Spiner::DataBox<
       double, Spiner::RegularGrid1D<double>,
       Spiner::InterpolationStencils::Rational>;

   // Host-side example: four samples of f(x) = x*x.
   RationalBox db(4);
   db.setRange(0, 0.0, 3.0, 4);
   for (int i = 0; i < 4; ++i) {
     db(i) = double(i) * double(i);
   }
   auto [value, df_dx] = db.interpToRealWithGrads(1.5);
   // value = 2.25, df_dx = 3.0
   db.finalize();

The policy supports fully interpolated rank-1 and rank-2 databoxes and
requires at least four points on every interpolated axis. All four built-in
grid types are supported. For 2D, call
``db.interpToRealWithGrads(y, x)``; the result contains ``value``, ``d_dx2``
(the y derivative), and ``d_dx1`` (the x derivative), in that order.

The one-dimensional stencil uses rational interpolation in interior
intervals and quadratic formulas in the first and last intervals. Strictly
outside the table, it uses linear secant extrapolation from the nearest
two endpoint samples. At an exact table endpoint, it returns the interior
derivative, which can differ from the exterior secant slope; derivatives
can therefore jump at the endpoints. These axis rules also apply in 2D.

The policy affects only ``interpToRealWithGrads``. Value-only
``interpToReal`` calls remain linear, even on a ``RationalBox``. To obtain
only the rationally interpolated value, use
``db.interpToRealWithGrads(x).value`` (or the two-coordinate overload).

Mixed interpolation and indexing
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In the case where an array has some dimensions that are discrete and
some that are interpolatable, one can fuse interpolation and indexing
into a single operation. These operations are still named
``DataBox::interpToReal``, but one of the input arguments is an
integer instead of a floating point number. The location of the
integer in the function signature indicates which dimension in the
``DataBox`` is indexed. For example:

.. cpp:function:: T DataBox::interpToReal(const T x3, const T x2, const T x1, const int idx) const;

interpolates the three slower-moving indices and indexes the fastest
moving index. On the other hand,

.. cpp:function:: T DataBox::interpToReal(const T x4, const T x3, const T x2, const int idx, const T x1) const;

interpolates the fastest moving index, then indexes the
second-fastest, then interpolates the remaining three slower. The
above fused operations are the only ones currently supported.

Interpolating into another ``DataBox``
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

There is limited functionality for filling a ``DataBox`` with the
interpolated values of another ``DataBox``. For example, the method

.. cpp:function:: void DataBox::interpFromDB(const DataBox &src, const T x);

fills the ``DataBox`` with the faster moving elements of ``src``
interpolated to ``x`` in the slowest-moving direction. Similarly,

.. cpp:function:: void DataBox::interpFromDB(const DataBox &src, const T x2, const T x1);

interpolates the two slowest-moving directions of ``src`` to ``x2``
and ``x1``.

.. warning::

  ``interpFromDB`` never allocates, reshapes, or frees memory. The
  ``DataBox`` being filled must already exist with storage and the shape
  of the remaining faster-moving dimensions of ``src``. In other words,
  it must have rank ``src.rank() - 1`` (or ``src.rank() - 2``) and
  ``dim(i) == src.dim(i)`` for each remaining dimension. Both ``DataBox``
  objects must be accessible from the execution space in which
  ``interpFromDB`` is called, for example both on device inside a
  kernel. 

Only values are written. The grids and index types of the ``DataBox``
being filled are left unchanged, so call ``setRange`` on it yourself if
you want to interpolate it later. For example:

.. code-block:: cpp

  // src has shape (NZ, NY, NX)
  Spiner::DataBox<double> dst(NY, NX);
  dst.interpFromDB(src, z);
  dst.setRange(1, src.range(1));
  dst.setRange(0, src.range(0));
  double val = dst.interpToReal(y, x);
  // ...
  free(dst);

.. warning::

  Be careful when the grids own their own memory, as
  ``NonUniformGrid1D`` does. ``setRange`` copies a grid's ownership
  along with it, so passing ``src.range(i)`` would make both ``DataBox``
  objects owners of the same grid memory. To share the grids, pass a
  non-owning handle instead:

  .. code-block:: cpp

    Spiner::NonUniformGrid1D<double> x_grid;
    x_grid.shallowCopy(src.range(0));
    dst.setRange(0, x_grid);

  Then only ``src`` frees the grid memory, through ``finalizeGrids``.

File I/O
^^^^^^^^^

If `hdf5`_ is enabled, ``Spiner`` can save an array to or load an
array from disk. Each array so-saved is also saved with the
``IndexType`` and independent variable ranges bundled with it, so that
knowledge of how to interpolate the data is automatically
available.

.. _`hdf5`: https://www.hdfgroup.org/solutions/hdf5/

The following methods are supported:

.. cpp:function:: herr_t DataBox::saveHDF(const std::string &filename) const;
   
saves the ``DataBox`` to a file with ``filename``.

.. cpp:function:: herr_t DataBox::saveHDF(hid_t loc, const std::string &groupname) const;

saves the ``DataBox`` as an hdf5 group at the location ``loc`` in an hdf5 file.

.. cpp:function:: DataBox::loadHDF(const std::string &filename);

fills the ``DataBox`` from information in the root of a file with ``filename``.

.. cpp:function:: DataBox::loadHDF(hid_t loc, const std::string &groupname);

fills the ``DataBox`` from information in the group with ``groupname``
based at location ``loc`` in the file.

.. warning::
  HDF5 I/O is only supported for single- and double-precision types at this time.

Miscellany
^^^^^^^^^^^

Here we list a few convenience functions available that were not
covered elsewhere.

.. cpp:function:: T DataBox::min() const;

and

.. cpp:function:: T DataBox::max() const;

compute and return the minimum and maximum values (respectively) in the array.

.. cpp:function:: int rank() const;

returns the rank (number of dimensions) of the array.

.. cpp:function:: int size() const;

returns the total number of elements in the underlying array.

.. cpp:function:: int sizeBytes() const;

returns the total size of the underlying array in bytes.

.. cpp:function:: int dim(int i) const;

returns the size in a given dimension/direction, indexed from one.

Generative AI was used to assist with modifications to this page.