# Introduction

Unidata is in the process of developing a Zarr [] based variant of netcdf. As part of this effort, it was necessary to implement some support for chunking. Specifically, the problem to be solved was that of extracting a hyperslab of data from an n-dimensional variable (array in Zarr parlance) that has been divided into chunks (in the HDF5 sense). Each chunk is stored independently in the data storage -- Amazon S3, for example.

The algorithm takes a series of R slices of the form (first,stop,stride), where R is the rank of the variable. Note that a slice of the form (first, count, stride), as used by netcdf, is equivalent because stop = first + count*stride. These slices form a hyperslab.

The goal is to compute the set of chunks that intersect the hyperslab and to then extract the relevant data from that set of chunks to produce the hyperslab.

It appears from web searches that this algorithm is nowhere documented in the form of a high level pseudo code algorithm. It appears to only exist in the form of code in HDF5, Zarr, and probably TileDB, and maybe elsewhere.

What follows is an attempt to reverse engineer the algorithm used by the Zarr code to compute this intersection. This is intended to be the bare-bones algorithm with no optimizations included. Thus, its performance is probably not the best, but it should work in all cases. Some understanding of how HDF5 chunking is used is probably essential for understanding this algorithm.

The original python code relies heavily on the use of Python iterators, which was might be considered a mistake. Iterators are generally most useful in two situations: (1) it makes the code clearer -- arguably false in this case, and (2) there is a reasonable probability that the iterators will be terminated before they end, thus improving efficiency. This is demonstrably false for the Zarr code.

This code instead uses the concept of an *odometer*,
which is a way to iterate over all the elements of an
n-dimensional object. Odometer code already exists in several
places in the existing netcdf-c library, among them
*ncdump/nciter.c*, *ncgen/odom.c*, and *libdap4/d4odom.c*. It
is also equivalent to the Python *itertools.product* iterator
function.

This algorithm assumes the following items of input:
1. variable (aka array) - multidimensional array of data values
2. dimensions - the vector of dimensions defining the rank of the array;
this set comes from the dimensions associated with the array;
so given v(d1,d2) where (in netcdf terms d1 and d2 are the
dimension names and where for example, d1=10, d2=20
the dimensions given to this algorithm are the ordered vector (d1,d2),
or equivalently (10,20).
3. Chunk sizes - for each dimension in the dimension set, there is defined
a chunk size along that dimension. Thus, there is also an ordered vector
of chunk sizes, (c1,c2) corresponding to the dimension vector (d1,d2).
4. slices - a vector of slice instances defining the subset of data to be
extracted from the array. A single slice as used here is of the form
(start,stop,stride) where start is the starting position with respect to
a dimension, stop is the last position + 1 to extract, and stride is
the number of positions to skip when extracting data. Note that this
is different than the netcdf-c nc*get*vars slices, which are of the form
(start,count,stride). The two are equivalent since stop = start+(count*stride).
When extracting data from a variable of rank R, one needs to specify R slices
where each slice corresponds to a dimension of the variable.

At a high-level, the algorithm works by separately analyzing
each dimension of the array using the corresponding slice and
corresponding chunk size. The result is a set of *projections*
specific to each dimension. By taking the cross-product of these
projections, one gets a vector of projection-vectors that can be
evaluated to extract a subset of the desired data for storage in
an output array.

It is important to note that this algorithm operates in two phases. In the first phase, it constructs the projections for each dimension. In the second phase, these projections are combined as a cross-product to provide subsetting for the true chunks, which are R-dimensional rectangles.

What follows is the algorithm is written as a set of pseudo-code procedures to produce the final output given the above inputs.

## Notations:

- floordiv(x,y) = floor((x / y))
ceildiv(x,y) = ceil((x / y))

## Notes:

The ith slice is matched to the ith dimension for the variable

- The ith slice is matched to the ith chunksize for the variable
- The zarr code uses iterators, but this code converts to using vectors and odometers for (one hopes) some clarity and consistency with existing netcdf-c code.

## Global Type Declarations

```
class Slice {
int start
int stop
int stride
}
class SliceIndex { // taken from zarr code see SliceDimIndexer
int chunk0 // index of the first chunk touched by the slice
int nchunks // number of chunks touched by this slice index
int count //total number of output items defined by this slice index
Projection projections[nchunks]; // There are multiple projections
// derived from the original slice:
// one for each chunk touched by
// the original slice
}
class Projection {
int chunkindex;
Slice slice; // slice specific to this chunk
int outpos; // start point in the output to store the extracted data
}
```

## Global variables

In order to keep argument lists short, certain values are assumed to be globally defined and accessible. * R - the rank of the variable * dimlen[R] - the length of the dimensions associated with the variable * chunklen[R] - the length of the chunk sizes associated with the variable * int zeros[R] - a vector of zeros * int ones[R] - a vector of ones

## Procedure EvaluateSlices

```
// Goal: Given the projections for each slice being applied to the
// variable, create and walk all possible combinations of projection
// vectors that can be evaluated to provide the output data
void EvaluateSlices(
Slice slices[R], // the complete set of slices
T output // the target storage for the extracted data: its type is T
)
{
int i;
SliceIndex allsliceindices[R];
Odometer odometer;
nchunks[R]; // the vector of nchunks from projections
int indices[R];
// Compute the slice index vector
allsliceindices = compute_all_slice_indices(slices);
// Extract the chunk0 and nchunks vectors
for(i=0;i<R;i++) {
nchunks[i] = allsliceindices[i].nchunks;
}
// Create an odometer to walk nchunk combinations
odometer = Odometer.new(R,zeros,nchunks); // iterate each "wheel[i]" over 0..nchunk[i] with R wheels
// iterate over the odometer: all combination of chunk indices in the projections
for(;odometer.more();odometer.next()) {
chunkindices = odometer.indices();
ApplyChunkIndices(chunkindices,output,allsliceindices);
}
}
```

## Procedure ApplyChunkIndices

```
// Goal: given a vector of chunk indices from projections,
// extract the corresponding data and store it into the
// output target
void ApplyChunkIndices(
int chunkindices[R], // indices chosen by the parent odometer
T output, // the target storage for the extracted data
SliceIndex allsliceindices[R]
)
{
int i;
SliceIndex subsliceindices[R];
chunk0[R]; // the vector of chunk0 values from projections
Projection projections[R];
int[R] start, stop,stride; // decomposed set of slices
int[R] outpos; // capture the outpos values across the projections
int ouputstart;
// This is complicated. We need to construct a vector of slices
// of size R where the ith slice is determined from a projection
// for the ith chunk index of chunkindices. We then iterate over
// that odometer to extract values and store them in the output.
for(i=0;i<R;i++) {
int chunkindex = chunkindices[i];
Slice slices[R];
projections[i] = allsliceindices[i].projections[chunkindex];
slices[i] = projections[i].slice;
outpos[i] = projections[i].outpos;
}
// Compute where the extracted data will go in the output vector
outputstart = computelinearoffset(R,outpos,dimlen);
GetData(slices,outputstart,output);
}
```

## Procedure GetData

```
// Goal: given a set of indices pointing to projections,
// extract the corresponding data and store it into the
// output target.
void GetData(
Slice slices[R],
int chunksize, // total # T instances in chunk
T chunk[chunksize],
int outputstart,
T output
)
{
int i;
Odometer sliceodom,
sliceodom = Odometer.new(R, slices);
// iterate over the odometer to get a point in the chunk space
for(;odom.more();odom.next()) {
int chunkpos = odometer.indices(); // index
}
}
```

## Procedure compute_all_slice_projections

```
// Goal:create a vector of SliceIndex instances: one for each slice in the top-level input
Projection[R]
compute_all_slice_projections(
Slice slice[R], // the complete set of slices
{
int i;
SliceIndex projections[R];
* for i in 0..R-1
* projections[i] = compute_perslice_projections(dimlen[i],chunklen[i],slice[i])
* return projections
```

## Procedure compute_perslice_projections

### Goal:

For each slice, compute a set of projections from it wrt a dimension and a chunk size associated with that dimension.

### Inputs:

dimlen -- dimension length

- chunklen -- chunk length associated with the input dimension
slice=(start,stop,stride) -- associated with the dimension

### Outputs:

Instance of SliceIndexs

### Data Structure:

````

### Computations:

- count = max(0, ceildiv((stop - start), stride))
- total number of output items defined by this slice (equivalent to count as used by nc
*get*vars)

- total number of output items defined by this slice (equivalent to count as used by nc
- nchunks = ceildiv(dim
*len, dim*chunk_len)- number of chunks touched by this slice

- chunk0 = floordiv(start,chunklen)
- index (in 0..nchunks-1) of the first chunk touched by the slice

- chunkn = ceildiv(stop,chunklen)
- index (in 0..nchunks-1) of the last chunk touched by the slice

- n = ((chunkn - chunk0) + 1)
- total number of touched chunks
- the index i will range over 0..(n-1)
- is this value the same as nchunks?

- For each touched chunk index we compute a projection specific to that chunk, hence there are n of them.
- projections.index[i] = i
- projections.offset[i] = chunk0 * i
- remember: offset is only WRT this dimension, not global

- projections.limit[i] = min(dimlen, (i + 1) * chunklen)
- end of this chunk but no greater than dimlen

- projections.len[i] = projections.limit[i] - projections.offset[i]
- actual limit of the ith touched chunk; should be same as chunklen except for last length because of the min function in computing limit[i]

- projections.start[i]:
- This is somewhat complex because for the first projection, the start is the slice start, but after that, we have to take into account that for a non-one stride, the start point in a projection may be offset by some value in the range of 0..(stride-1)
- i == 0 => projections.start[i] = start - projections.offset[i]
- initial case the original slice start is within the first projection

- i > 0 => projections.start[i] = start - projections.offset[i]
- prevunused[i] = (projections.offset[i] - start) % stride
- prevunused[i] is an intermediate computation and need not be saved
- amount unused in previous chunk => we need to skip (stride-prevunused[i]) in this chunk

- prevunused[i] > 0 => projections.start[i] = stride - prevunused[i]

- prevunused[i] = (projections.offset[i] - start) % stride

- projections.stop[i]:
- stop > projections.limit[i] => projections.stop[i] = projections.len[i]
- stop <= projections.limit[i] => projections.stop[i] = stop - projections.offset[i]
- selection ends within current chunk

- projections.outpos[i] = ceildiv(offset - start, stride)
- "location" in the output array to start storing items; again, per slice, not global

## Procedure computelinearoffset(outpos,dimlen);

```` // Goal: Given a set of per-dimension indices, compute the corresponding linear position.

```
int
computelinearoffset(int R,
int outpos[R],
int dimlen[R]
)
{
int offset;
int i;
offset = 0;
for(i=0;i<R;i++) {
offset *= dimlen[i];
offset += outpos[i];
}
return offset;
}
```

# Appendix: Odometer Code

```
class Odometer
{
int R; // rank
int start[R];
int stop[R]
int stride[R];
int index[R]; // current value of the odometer
procedure new(int R, int start[R], int stop[R]) { return new(R, start,stop,ones);}
procedure new(int rank, Slice slices[R])
{
int i;
int start0[R];
int stop0[R];
int stride0[R];
R = rank;
for(i=0;i<R;i++) {
start = slices[i].start;
stop = slices[i].stop;
stride = slices[i].stride;
}
for(i=0;i<R;i++) {index[i] = start[i];}
}
boolean
procedure more(void)
{
return (index[0] < stop[0]);
}
procedure next(void)
{
int i;
for(i=R-1;i>=0;i--) {
index[i] += stride[i];
if(index[i] < stop[i]) break;
if(i == 0) break; // leave the 0th entry if it overflows
index[i] = start[i]; // reset this position
}
}
// Get the value of the odometer
int[R]
procedure indices(void)
{
return indices;
}
}
```