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
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.
This is the approach taken by the
Those axes are represented by the
As usual we construct the
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
Here we first check that the axis is a nonempty 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.
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
Finally, we can copy another gridlike object from its
Converting an
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
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
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.
Program 2 uses this approach to print successive states of an array of three indices as they are iterated from zero to three.
We could generalise this by passing in a function to call with the array of indices as its argument, much like the
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.
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.
The
Program 5 uses this function to iterate over the nodes of an
Whilst we're about it we might as well implement a way to iterate over multiple indices in reverse, as done by
Finally, the implementation of the
[2] Cubic Line Division, www.thusspakeak.com, 2018.
[3] We're Not For Turning, www.thusspakeak.com, 2018.
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 ofak.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 : lr;}); for(i=1;i<n;++i) { if(axis[i]===axis[i1]) { throw new Error('duplicate axis value in ak.grid'); } } }
Here we first check that the axis is a nonempty 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 gridlike 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 looplike 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[n1]<u[n1])) {n; a[n]=l[n];} return n>0; }; nextState[ak.ARRAY_T][ak.NUMBER_T] = function(n, a, l, u) { while(n>0 && !(++a[n1]<u)) {n; a[n]=l[n];} return n>0; }; nextState[ak.ARRAY_T][ak.UNDEFINED_T] = function(n, a, u) { while(n>0 && !(++a[n1]<u[n1])) {n; a[n]=0;} return n>0; }; nextState[ak.NUMBER_T][ak.ARRAY_T] = function(n, a, l, u) { while(n>0 && !(++a[n1]<u[n1])) {n; a[n]=l;} return n>0; }; nextState[ak.NUMBER_T][ak.NUMBER_T] = function(n, a, l, u) { while(n>0 && !(++a[n1]<u)) {n; a[n]=l;} return n>0; }; nextState[ak.NUMBER_T][ak.UNDEFINED_T] = function(n, a, u) { while(n>0 && !(++a[n1]<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[n1]>=l[n1])) {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[n1]>=l[n1])) {n; a[n]=u1;} return n>0; }; prevState[ak.ARRAY_T][ak.UNDEFINED_T] = function(n, a, u) { while(n>0 && !(a[n1]>=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[n1]>=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[n1]>=l)) {n; a[n]=u1;} return n>0; }; prevState[ak.NUMBER_T][ak.UNDEFINED_T] = function(n, a, u) { while(n>0 && !(a[n1]>=0)) {n; a[n]=u1;} 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