Cuboid Space Division

| 0 Comments

Over the last few months we have been taking a look at algorithms for interpolating over a set of points \(\left(x_i,y_i\right)\) in order to approximate values of \(y\) between the nodes \(x_i\). We began with linear interpolation[1] which connects the points with straight lines and is perhaps the simplest interpolation algorithm. Then we moved on to cubic spline interpolation[2] which yields a smooth curve by specifying gradients at the nodes and fitting cubic polynomials between them that match both their values and their gradients. Next we saw how this could result in curves that change from increasing to decreasing, or vice versa, between the nodes and how we could fix this problem by adjusting those gradients[3].
I concluded by noting that, even with this improvement, the shape of a cubic spline interpolation is governed by choices that are not uniquely determined by the points themselves and that linear interpolation is consequently a more mathematically appropriate scheme, which is why I chose to generalise it to other arithmetic types for \(y\), like complex numbers or matrices, but not to similarly generalise cubic spline interpolation.

The obvious next question is whether or not we can also generalise the nodes to other arithmetic types; in particular to vectors so that we can interpolate between nodes in more than one dimension.
The simplest way to do this is to require that the nodes sit at the intersection of marks upon a set of axes rather than at those of a single one. For example, in two dimensions that they are vectors \(\mathbf{x}_{ij}\) defined as
\[ \mathbf{x}_{ij} = \left(x_{0i}, x_{1j}\right) \]
for some sets of scalars \(x_{0i}\) and \(x_{1j}\) so that they lie at the corners of a set of rectangles that cover the area bounded by the smallest and largest positions of the nodes along both axes. In three dimensions, that they lie at the vertices of boxes, or rectangular cuboids, that cover the volume bounded by the extreme positions on all three axes and, in more dimensions, at the vertices of their hyperdimensional analogues.

Now, we'll need a way to represent such nodes before we can even begin to think about implementing an algorithm to interpolate between them and so we might as well get on with that now and worry about the algorithmic details later.

Boxing Clever

If we were to represent the nodes with an array of ak.vector objects then we would be faced with the tedious business of checking that they occupy every intersection of the marks upon the axes. A much better scheme is to keep track of the marks upon the axes and create ak.vector nodes on the fly so that they're in the correct positions by construction.

This is the approach taken by the ak.grid type given in listing 1, which is essentially a container of marked axes.

Listing 1: ak.grid
ak.GRID_T = 'ak.grid';

function Grid(){}
Grid.prototype = {TYPE: ak.GRID_T, valueOf: function(){return ak.NaN;}};

ak.grid = function() {
 var g = new Grid();
 var state = [];
 var arg0  = arguments[0];

 constructors[ak.nativeType(arg0)](state, arg0, arguments);

 g.dims = function() {return state.length;};
 g.length = function(i) {return state[i].length;};
 g.lengths = function() {return state.map(function(a){return a.length;});};
 g.axis = function(i) {return state[i].slice(0);};
 g.axes = function() {return state.map(function(a){return a.slice(0);});};
 g.at = function(i, j) {return state[i][j];};
 g.map = function(a) {return mapAxes(state, a);};

 g.toArray  = g.axes;
 g.toString = function() {return toString(state);};

 g.toExponential = function(d) {
  return toString(state, function(x){return x.toExponential(d);});
 };

 g.toFixed = function(d) {
  return toString(state, function(x){return x.toFixed(d);});
 };

 g.toPrecision = function(d) {
  return toString(state, function(x){return x.toPrecision(d);});
 };

 return Object.freeze(g);
};

var constructors = {};

Those axes are represented by the state array which, after construction, will contain arrays representing the tick marks upon each axis, as is strongly hinted at by the basic property access methods. Specifically dims returns the number of axes, or equivalently the dimension of the space that they define, length returns the number of tick marks on the \(i^\mathrm{th}\) axis and lengths an array of the number of them on each axis, axis returns a copy of the \(i^\mathrm{th}\) axis and axes an array containing copies of all of them, and finally at returns the \(j^\mathrm{th}\) tick mark on the \(i^\mathrm{th}\) axis.

As usual we construct the ak.grid by dispatching to a constructors object, with the first target creating a set of identical axes from the first argument specifying their number and the second their ticks, as shown by listing 2.

Listing 2: Identical Axes
constructors[ak.NUMBER_T] = function(state, n, args) {
 var axis = args[1];

 if(n<=0 || n!==ak.floor(n)) {
  throw new Error('invalid dimension in ak.grid');
 }
 if(ak.nativeType(axis)!==ak.ARRAY_T) {
  throw new Error('invalid axis in ak.grid');
 }
 axis = axis.slice(0);
 sortAxis(axis);

 state.length = n;
 while(n-->0) state[n] = axis;
};

Note that we don't bother dispatching on the second argument since we require it to be an array and so instead throw an exception if it isn't one. Aside from this we simply check that the number of axes is a positive integer and set the elements of an array of that length to a copy of the axis that has been processed by the sortAxis function defined in listing 3.

Listing 3: Sorting The Axes
function sortAxis(axis) {
 var n = axis.length;
 var i;

 if(n===0) throw new Error('empty axis in ak.grid');

 for(i=0;i<n;++i) {
  axis[i] = Number(axis[i]);
  if(isNaN(axis[i])) {
   throw new Error('invalid axis value in ak.grid');
  }
 }
 axis.sort(function(l,r){return l===r ? 0 : l-r;});
 for(i=1;i<n;++i) {
  if(axis[i]===axis[i-1]) {
   throw new Error('duplicate axis value in ak.grid');
  }
 }
}

Here we first check that the axis is a non-empty array of numbers, or at least of values that can be converted to numbers, before sorting them in ascending order, taking care to correctly handle infinite values by explicitly returning zero if both arguments of the comparison are equal. Finally, we ensure that the tick marks are unique by throwing an exception if any adjacent elements of the sorted array are equal.

To construct a grid with different axes we pass an array of arrays to the constructor which dispatches to the function defined in listing 4.

Listing 4: Different Axes
constructors[ak.ARRAY_T] = function(state, axes) {
 var n = axes.length;
 if(n===0) throw new Error('empty axes in ak.grid');

 state.length = n;
 while(n-->0) {
  if(ak.nativeType(axes[n])!==ak.ARRAY_T) {
   throw new Error('invalid axis in ak.grid');
  }
  state[n] = axes[n].slice(0);
  sortAxis(state[n]);
 }
};

The only difference between this and the previous constructor is that we must validate, sort and copy the axis one at a time from the argument to the state.

Finally, we can copy another grid-like object from its dims, length and at property access methods using the constructor given in listing 5, also accepting a numeric property for dims.

Listing 5: Copying A Grid Object
constructors[ak.OBJECT_T] = function(state, grid) {
 var n = grid.dims;
 var m, axis;

 n = (ak.nativeType(n)===ak.FUNCTION_T) ? Number(n()) : Number(n);
 if(n<=0 || n!==ak.floor(n)) {
  throw new Error('invalid dimension in ak.grid');
 }

 state.length = n;
 while(n-->0) {
  m = Number(grid.length(n));
  if(m<=0 || m!==ak.floor(m)) {
   throw new Error('invalid axis length in ak.grid');
  }
  axis = new Array(m);
  while(m-->0) axis[m] = grid.at(n, m);
  sortAxis(axis);
  state[n] = axis;
 }
};

Converting an ak.grid to an array is trivially identical to copying its axes and so the toArray conversion method is simply a synonym for the axes property access method. Converting it to a string is rather more convoluted and so we forward to the toString helper function given in listing 6 from toString, toExponential, toFixed and toPrecision.

Listing 6: Converting A Grid To A String
function toString(axes, f) {
 var s = [];
 var n = axes.length;
 var i, j, a, m;

 s.push('[');

 for(i=0;i<n;++i) {
  a = axes[i];
  m = a.length;

  s.push('[');

  for(j=0;j<m;++j) {
   s.push(f ? f(a[j]) : a[j]);
   s.push(',');
  }

  s.pop();
  s.push(']');
  s.push(',');
 }

 s.pop();
 s.push(']');

 return s.join('');
}

This represents the grid as a square bracketed, comma separated list of axes which are represented as square bracketed, comma separated lists of the positions of the tick marks, removing the trailing comma before adding each closing bracket.

Finally, we come to the on the fly creation of the nodes with the map method which forwards the state and its indexing a argument to the mapAxes function given in listing 7.

Listing 7: Constructing A Node
function mapAxes(axes, index) {
 var n = axes.length;
 if(ak.nativeType(index)!==ak.ARRAY_T || index.length!==n) {
  throw new Error('invalid index in ak.grid.map');
 }
 index = index.slice(0);
 while(n-->0) index[n] = axes[n][index[n]];
 return ak.vector(index);
}

This assumes that the index is an array of the indices of the tick marks upon each axis at which the node sits and simply replaces them with their respective positions in a copy of the array, constructing an ak.vector with the result to represent the node. The mapping of array indices to vector nodes is demonstrated by program 1.

Program 1: Mapping The Nodes

Using nested loops to iterate over the nodes isn't exactly convenient since we'd need to rewrite the code for grids with different numbers of dimensions and so it's worth taking some time to develop a mechanism that enables us to avoid doing so before we move on to the business of interpolation.

Recursion To The Rescue!

The typical way to deal with such problems is to use recursion; specifically to create a function that iterates over an element of an array of indices calling itself at each step indicating the next element to iterate over, or takes some action with the values of the indices if there aren't any more to iterate over.
Program 2 uses this approach to print successive states of an array of three indices as they are iterated from zero to three.

Program 2: Recursing Over Indices

We could generalise this by passing in a function to call with the array of indices as its argument, much like the forEach method of JavaScript's Array type, but I'd prefer to separate the process of iteration from the action to take at each step, not least because it's trivial to implement forEach-like behaviour from loop-like behaviour but not so much the other way around.
To facilitate looping we need a function that updates the array of indices at each step and returns a boolean to indicate when the iteration should end. We can do this by incrementing the last of them, returning true if it's less than the upper limit, setting it to zero and recursively incrementing the next last index if it's not and returning false if the first index has reached its limit, as demonstrated by program 3.

Program 3: Loopy Recursion

Now we don't actually need to use recursion to update the array of indices; instead we can simply increment the last, return true if it's less than the upper bound, set it to zero and iterate to the previous index if not, returning false if we've passed the first index, as illustrated by program 4.

Program 4: Loopy Iteration

The ak.nextState function given in listing 8 uses this scheme to iterate over an array of indices between single bounds, passed as numbers, and multiple bounds, passed as arrays, distinguishing between them by dispatching to functions based upon the types of the arguments and defaulting to lower bounds of zero.

Listing 8: ak.nextState
var nextState = {};
nextState[ak.ARRAY_T] = {};
nextState[ak.NUMBER_T] = {};

nextState[ak.ARRAY_T][ak.ARRAY_T] = function(n, a, l, u) {
 while(n>0 && !(++a[n-1]<u[n-1])) {--n; a[n]=l[n];}
 return n>0;
};

nextState[ak.ARRAY_T][ak.NUMBER_T] = function(n, a, l, u) {
 while(n>0 && !(++a[n-1]<u)) {--n; a[n]=l[n];}
 return n>0;
};

nextState[ak.ARRAY_T][ak.UNDEFINED_T] = function(n, a, u) {
 while(n>0 && !(++a[n-1]<u[n-1])) {--n; a[n]=0;}
 return n>0;
};

nextState[ak.NUMBER_T][ak.ARRAY_T] = function(n, a, l, u) {
 while(n>0 && !(++a[n-1]<u[n-1])) {--n; a[n]=l;}
 return n>0;
};

nextState[ak.NUMBER_T][ak.NUMBER_T] = function(n, a, l, u) {
 while(n>0 && !(++a[n-1]<u)) {--n; a[n]=l;}
 return n>0;
};

nextState[ak.NUMBER_T][ak.UNDEFINED_T] = function(n, a, u) {
 while(n>0 && !(++a[n-1]<u)) {--n; a[n]=0;}
 return n>0;
};

ak.nextState = function(a, l, u) {
 return nextState[ak.nativeType(l)][ak.nativeType(u)](a.length, a, l, u);
};

Program 5 uses this function to iterate over the nodes of an ak.grid object by using its lengths property as the upper bound.

Program 5: Iterating Over The Nodes

Whilst we're about it we might as well implement a way to iterate over multiple indices in reverse, as done by ak.prevState in listing 9.

Listing 9: ak.prevState
var prevState = {};
prevState[ak.ARRAY_T] = {};
prevState[ak.NUMBER_T] = {};

prevState[ak.ARRAY_T][ak.ARRAY_T] = function(n, a, l, u) {
 while(n>0 && !(--a[n-1]>=l[n-1])) {--n; a[n]=u[n]-1;}
 return n>0;
};

prevState[ak.ARRAY_T][ak.NUMBER_T] = function(n, a, l, u) {
 while(n>0 && !(--a[n-1]>=l[n-1])) {--n; a[n]=u-1;}
 return n>0;
};

prevState[ak.ARRAY_T][ak.UNDEFINED_T] = function(n, a, u) {
 while(n>0 && !(--a[n-1]>=0)) {--n; a[n]=u[n]-1;}
 return n>0;
};

prevState[ak.NUMBER_T][ak.ARRAY_T] = function(n, a, l, u) {
 while(n>0 && !(--a[n-1]>=l)) {--n; a[n]=u[n]-1;}
 return n>0;
};

prevState[ak.NUMBER_T][ak.NUMBER_T] = function(n, a, l, u) {
 while(n>0 && !(--a[n-1]>=l)) {--n; a[n]=u-1;}
 return n>0;
};

prevState[ak.NUMBER_T][ak.UNDEFINED_T] = function(n, a, u) {
 while(n>0 && !(--a[n-1]>=0)) {--n; a[n]=u-1;}
 return n>0;
};

ak.prevState = function(a, l, u) {
 return prevState[ak.nativeType(l)][ak.nativeType(u)](a.length, a, l, u);
};

Finally, the implementation of the ak.grid type and the ak.nextState and ak.prevState functions can be found in Grid.js, NextState.js and PrevState.js respectively.

References

[1] Chalk The Lines, www.thusspakeak.com, 2018.

[2] Cubic Line Division, www.thusspakeak.com, 2018.

[3] We're Not For Turning, www.thusspakeak.com, 2018.

Leave a comment