Growing a Tabulation¶

Many quantities in Galacticus are expensive enough that they are tabulated once and interpolated thereafter. The range which must be tabulated is usually not known in advance: it is discovered from the values actually requested as a model runs. A tabulation therefore has to grow, and the question this page answers is how to grow one without changing the values it has already returned.

The mechanism is an absolute lattice: a set of abscissae fixed once and for all by the gridding scheme, entirely independent of which range any particular table happens to span. A table built on such a lattice can be extended, and every point it already held—and every value interpolated between those points—is unchanged, bit for bit. The pieces are Range_Pinned and rangeLattice (in Numerical_Ranges), the extend methods of the table1D and table2DLogLogLin classes, and, for tabulations expensive enough to keep between runs, Table_Caches.

Why a naively-grown tabulation is a problem¶

The natural way to grow a table is to widen its bounds until they enclose the new request, redistribute the requested number of points over the widened range, and paste the previously computed values back in. This does not work, and the reason is worth understanding before reaching for the machinery which avoids it.

The bounds so obtained depend on the history of the requests. They are grown by repeated multiplication by some factor from an origin set by whichever value happened to be asked for first, so two runs which end up covering the same range of values, but reach it in a different order, tabulate at different abscissae. Nothing about the model has changed, yet every interpolated value differs.

Worse, the abscissae pasted back in no longer belong to the table which now holds them. Redistributing \(N\) points over the widened range gives a new spacing, but the preserved values retain the abscissae they were computed at, so the table becomes a patchwork of spacings which no single interpolation factor describes. Where such a table is indexed arithmetically—from a stored inverse spacing, rather than by searching—the index computed drifts further from the truth with each extension. This is not hypothetical: it is exactly the failure which the expansion-factor tabulation in source/cosmology/functions/matter_lambda.F90 used to suffer, where after sufficiently many extensions the computed index fell beyond the end of the table.

The absolute lattice¶

A lattice is specified by a gridding scheme and a density of points, pointsPer \(=N\). Its points are the \(x_j\), for all integers \(j\):

gridSchemePerDecade — \(x_j = 10^{j/N}\);
gridSchemePerOctave — \(x_j = 2^{j/N}\);
gridSchemePerUnit — \(x_j = j/N\).

Since \(x_j\) is a function of \(j\) and \(N\) alone, the lattice does not know, and cannot depend on, the range over which any particular table is built. Two tables built on the same lattice therefore have identical abscissae wherever their ranges overlap.

A rangeLattice object names a contiguous stretch of such a lattice: the count points with integer indices indexMinimum through indexMinimum+count-1. It carries no abscissae of its own—they are regenerated from the integer indices on demand by values() (or valuesLogarithmic()), so they are bit-identical to those of any other range on the same lattice. Its useful methods are:

isDefined() — is this a usable lattice? A default-initialized rangeLattice is not, and that is how “nothing is tabulated yet” is represented;
covers(lattice) — is lattice on the same lattice as this one, with every one of its points also a point of this one? This is the test for “is the existing tabulation already sufficient?”;
minimum(), maximum(), indexMaximum(), values(), valuesLogarithmic();
step() (perUnit only) and stepLogarithmic() (logarithmic schemes only) — the spacing of the points.

step() and stepLogarithmic() are evaluated as pure functions of pointsPer, and deliberately not as the difference of two neighboring lattice points: that difference varies in its final bits with position along the lattice, so a table which derived its interpolation factor from it would change every interpolated value by of order one unit in the last place whenever its lower bound moved—defeating the very reuse the lattice exists to provide. Always take the spacing from the lattice.

Pinning a range to the lattice¶

Range_Pinned turns a request—one value, or an array of values, which the tabulation must cover—into the rangeLattice to tabulate:

use :: Numerical_Ranges, only : Range_Pinned, rangeLattice, gridSchemePerDecade
...
type(rangeLattice) :: lattice
...
lattice=Range_Pinned(                                    &
     &                              time               , &
     &                              pointsPerDecade    , &
     &                              gridSchemePerDecade, &
     &               marginFactor  =2.0d0              , &
     &               anchorEvery   =pointsPerDecade/2  , &
     &               latticeCurrent=self%lattice         &
     &              )

The bounds are found by, in order:

applying the safety margin to the range of the requested values: a multiplicative marginFactor (default \(2\)) at each end for the logarithmic schemes, or an additive marginOffset (default \(1\)) for perUnit. The margin means that a request which creeps outward by a little does not force a retabulation for every step it takes;
taking the union with rangeCurrent, if given—a pair of values which the range must span, used where some bound must always be included but is not itself known to lie on the lattice;
pinning the bounds outward to lattice points whose index is a multiple of anchorEvery;
clamping to the hard limits limitMinimum and limitMaximum, if given, with the clamped edge snapped inward to the lattice so that all points remain lattice points;
taking the union with latticeCurrent, if given—in exact integer arithmetic, and after the clamping, so that a lattice already in use is never shrunk and extension of a table built on it can never fail.

`anchorEvery`¶

anchorEvery is the interval, in lattice steps, to which the bounds are pinned. It controls only how far beyond the request the table extends, and is independent of pointsPer, which controls the accuracy of interpolation within it. It defaults to pointsPer—whole decades, octaves, or unit intervals—which gives the coarsest set of possible ranges and hence the strongest determinism. A value of \(1\) pins to the lattice points themselves.

Intermediate values are for cases where a whole decade of extra tabulation is extravagant. On a cosmic time axis a decade spans most of the history of the universe, so anchorEvery=pointsPer/2 (half decades, a factor of \(\sqrt{10}\)) is the usual choice there; where each additional point costs an ODE integration, anchorEvery=1 may be the only affordable choice. Determinism, exact reuse on extension, and the mergeability of cached tables hold for any choice—only the granularity of the set of possible ranges changes.

`rangeCurrent` versus `latticeCurrent`¶

These two look similar and are not interchangeable. The distinction is the most common source of error in using Range_Pinned, so it is worth stating plainly:

the request is the target, and only the target receives the safety margin;
the range already tabulated is taken in union through latticeCurrent, after the margin and after the hard limits.

Folding the current bounds into the target instead—taking a min/max against them before calling Range_Pinned—applies the safety margin to a bound to which the margin has already been applied once. The range then ratchets outward by another factor on every retabulation, growing without limit. latticeCurrent exists precisely so that the current range can be honored without the margin being applied to it a second time.

rangeCurrent is for the different case of a bound which must always be spanned but does not lie on the lattice: the present day, for a tabulation of integrals evaluated to the present day, or the current range of a table which was not itself built on a lattice. It is applied after the margin, so it raises the upper bound without also lowering the lower one.

Deciding whether to retabulate¶

Test the pinned lattice, not the requested value against the current bounds:

lattice=Range_Pinned(x,pointsPer,scheme,latticeCurrent=self%table%lattice)
if (self%table%lattice%covers(lattice)) return

Testing the request against the bounds would apply the safety margin only when the request happened to arrive outside them, so the range finally reached would depend on the order in which values were asked for—the very dependence pinning exists to remove. Pinning first and comparing lattices costs one cheap integer comparison and has no such history.

Extending the tabulation¶

extend grows a table onto a new lattice, creating it if it does not yet exist, and returns an isComputed mask which is true at those points whose values were carried over and false at those the caller must now evaluate:

logical, allocatable, dimension(:) :: isComputed
...
call self%table%extend(lattice,isComputed)
do i=1,lattice%count
   if (isComputed(i)) cycle
   call self%table%populate(expensiveFunction(self%table%x(i)),i)
end do

Two things to note. Take the abscissae from the table (or from lattice%values()), never recompute them independently, so that what is evaluated is exactly what will later be interpolated. And skipping the computed points is the entire purpose of the exercise—the mask is not advisory.

extend and x are methods of the abstract table1D, but populate is defined on the concrete table types. Where the table is held as a class(table1D), allocatable—as it must be to use Table_Caches—the loop therefore needs a select type around it:

select type (table_ => self%table)
type is (table1DLogarithmicLinear)
   do i=1,lattice%count
      if (isComputed(i)) cycle
      call table_%populate(expensiveFunction(table_%x(i)),i)
   end do
end select

The new lattice must be commensurate with (same scheme, same pointsPer) and must contain the lattice on which the table is currently tabulated; anything else is a fatal error. Range_Pinned guarantees this when it is passed the table’s current lattice as latticeCurrent.

extend is supported for the table types with uniformly-spaced internal abscissae—the linear and logarithmic types, with linear or cubic-spline interpolation—and is an error for the others. The two-dimensional table2DLogLogLin takes a lattice per axis, each pinned independently:

call self%table%extend(latticeX,latticeY,isComputed)

Here the preserved values form a rectangular block, so isComputed is true over that block and false over the L-shaped remainder. Both axes are extended in one call because the block preserved is the rectangle common to the two tabulations.

Tabulations held in raw arrays¶

Not every tabulation lives in a table object. For a plain rank-1 allocatable array, Range_Lattice_Extend does the whole job:

use :: Numerical_Ranges, only : Range_Lattice_Extend
...
call Range_Lattice_Extend(self%lattice,latticeNew,self%values,isComputed)
self%lattice=latticeNew

Range_Lattice_Extend_2D does the same for a rank-2 array, extending both axes at once. Where the tabulation is of some other rank, or the tabulated axis is not the first, use Range_Lattice_Offset: it returns the index offset at which the points of the current lattice appear within the new one, and leaves the caller to move its own arrays—typically a Move_Alloc followed by a single array-section assignment.

What is, and is not, preserved¶

The reuse guarantee is strong, and it is worth being precise about its scope.

Preserved bit-for-bit. The abscissae, which are regenerated from integer lattice indices rather than carried over, and the tabulated values, which are copied into the index window they occupy in the new table. That window is computed from the integer lattice indices, so no searching or floating-point comparison of abscissae is involved. For a linearly-interpolated table the interpolation is local and the spacing comes from the lattice, so every interpolated value on the overlap is unchanged too.

This matters to any consumer which integrates across a range of the tabulated variable, or compares quantities evaluated either side of a retabulation. Absent the guarantee, such a consumer sees the tabulation shift beneath it: a quantity computed before the table grew and one computed after are no longer strictly comparable, and differences which ought to be zero are not.

Not preserved: cubic-spline coefficients. A cubic spline is not local—every coefficient depends on every tabulated value—so adding points at either end changes the interpolant everywhere. The coefficients are therefore discarded and rebuilt. A cubic-spline table preserves its tabulated values exactly, but not the values it interpolates between them. Callers must accept this; where it is unacceptable, use a linearly-interpolated table.

A note on why Lattice_Value, which converts a lattice coordinate back to a value, is marked noinline: were it inlined into an array-filling loop the compiler would be free to vectorize the exponentiation, and a vectorized pow() need not agree in its final bits with the scalar one. The same lattice point could then differ between two tables merely because the loops filling them had different trip counts. Keeping it an opaque call forces every lattice point, everywhere, through the same scalar evaluation.

When carry-over is not valid¶

Carrying values over is correct only when each tabulated value depends on its own abscissa and on nothing else about the range. That holds for a function evaluated pointwise, and for an integral whose limits do not move—an integral from each tabulated time to the present day, for instance, is unaffected when the earliest time tabulated changes.

It does not hold when the values are generated by marching across the table. The expansion-factor tabulation in matter_lambda.F90 is the worked example: the expansion factor is integrated forward from an initial condition imposed at the earliest tabulated epoch, so every value in the table depends on where that integration began. Extend the table downward and the initial condition moves; every value which was carried over is then the solution of a different initial value problem from the one the new table represents.

The idiom for such a case is to discard everything by presenting an undefined lattice to the extension, which allocates afresh with isComputed false throughout:

carryOver=self%initialized .and. self%lattice%isDefined()
if (carryOver) carryOver=latticeNew%indexMinimum == self%lattice%indexMinimum
if (.not.carryOver) then
   if (allocated(self%values)) deallocate(self%values)
   self%lattice=rangeLattice()
end if
call Range_Lattice_Extend(self%lattice,latticeNew,self%values,isComputed)

Note that the test is on the integer lattice index of the lower bound, not on the value of the bound—exact, and free of any tolerance.

Where a tabulation is of this kind it is worth seeding its lower bound to a fixed value rather than letting it follow the request, so that the bound moves rarely, if ever, and the discard is rarely triggered.

Persistent caches¶

Where a tabulation is expensive enough to be worth keeping between runs, Table_Caches (source/objects/tables/caches.F90) writes it to a file under pathTypeDataDynamic. Because the table is built on an absolute lattice, a cached file whose range does not cover a new request need not be discarded: the cached abscissae are identical to the new ones wherever the two overlap, so the file can be merged with the table in memory, only the genuinely missing points computed, and the union written back. Cache files therefore grow monotonically toward the union of every range ever requested, rather than thrashing between runs.

The file name comes from Table_Cache_File_Name:

self%fileName   =Table_Cache_File_Name(                                                                                   &
     &                                 subDirectory    ='intergalacticMedium'                                           , &
     &                                 objectType      =char(self%objectType      (                                   )), &
     &                                 hashedDescriptor=char(self%hashedDescriptor(includeSourceDigest         =.true.  , &
     &                                                                             includeFileModificationTimes=.true.))  &
     &                                )
call Directory_Make(File_Path(self%fileName))

Passing includeSourceDigest=.true. and includeFileModificationTimes=.true. makes the name change whenever either the parameters of the object or the source code which generates the tabulation changes—so a stale cache is never read back. Build the name with this function rather than by hand: doing it by hand is how a discriminator comes to be omitted, and two different tabulations come to share a file.

The full idiom is then:

lattice=Range_Pinned(x,pointsPer,scheme,latticeCurrent=self%table%lattice)
if (self%table%lattice%covers(lattice)) return
! Merge in any tabulation already cached on disk, then re-evaluate the range required in the light of it.
call Table_Cache_Restore(self%table,self%fileName,status)
lattice=Range_Pinned(x,pointsPer,scheme,latticeCurrent=self%table%lattice)
call self%table%extend(lattice,isComputed)
if (any(.not.isComputed)) then
   ! Compute only the missing points. Note that no lock is held while doing so.
   select type (table_ => self%table)
   type is (table1DLogarithmicLinear)
      do i=1,lattice%count
         if (isComputed(i)) cycle
         call table_%populate(expensiveFunction(table_%x(i)),i)
      end do
   end select
   ! Merge with, and write back, the cache file.
   call Table_Cache_Store(self%table,self%fileName)
end if

Range_Pinned is called twice, and must be: restoring the cache may widen the table’s lattice, and the range finally tabulated must be the union of that with the request.

Locking is handled internally, and the important property is what is not locked. Table_Cache_Restore takes a shared lock, and Table_Cache_Store an exclusive one, but each holds it only for the duration of the file input/output—never across the tabulation itself, which is the expensive part. Table_Cache_Store re-reads the file under its exclusive lock before writing, so a range written by another process since our own read is merged rather than overwritten. Note that Table_Cache_Store may extend the table, since it returns holding the union of the cached and in-memory tabulations.

Merging in two dimensions is more restrictive than in one: the union of two rectangles is a rectangle only if one contains the other, so unless the cached tabulation covers that in memory (in which case it is adopted) or is covered by it (in which case there is nothing to do), the cached tabulation is reported and ignored.

Checklist¶

When adding or converting a tabulation:

Store a rangeLattice alongside the tabulation, or use the lattice component the table classes already provide.
Pin the range with Range_Pinned, passing the request as the target and the current lattice as latticeCurrent. Never fold the current bounds into the target.
Choose anchorEvery for the cost of over-tabulating, and pointsPer for the accuracy of interpolation. They are independent.
Return early when covers says the existing tabulation is sufficient.
Extend, and evaluate only the points for which isComputed is false.
Take abscissae, and any interpolation spacing, from the lattice or the table.
Check that carry-over is actually valid for the quantity being tabulated; if it is not, discard explicitly by resetting the lattice.
If the tabulation is expensive enough to keep between runs, wrap it in Table_Cache_Restore/Table_Cache_Store.

The unit tests in source/tests/make_ranges.F90 and source/tests/table_caches.F90 exercise all of this, and are a useful source of small, concrete examples.