We have spent a few months looking at how we might interpolate between sets of points \(\left(x_i, y_i\right)\), where the \(x_i\) are known as nodes and the \(y_i\) as values, to approximate values of \(y\) for values of \(x\) between the nodes, either by connecting them with straight lines^{[1]} or with cubic curves^{[2][3]}.
Last time^{[4]}, in preparation for interpolating between multidimensional vector nodes, we implemented the
With this in place we're ready to take a look at one of the simplest multidimensional interpolation schemes; multilinear interpolation.
The easiest way to describe this algorithm is to start with the special case of two dimensional nodes, known as bilinear interpolation, in which the nodes sit at the corners of a set of rectangles that divide up the area bounded by the four nodes that sit on the first and last tick marks on both axes.
To interpolate within one of those rectangles we can first choose a pair of opposing edges and draw a line perpendicular to them through the point of interest, then use linear interpolation along both of those edges to approximate the values at their points of intersection with that line and again between them to approximate the value at the point itself, as illustrated by figure 1 which marks the point of interest in red, the nodes in green and the points of intersection in black.
In three dimensions the nodes will sit at the corners of rectangular cuboids or, to put it more plainly, boxes. To interpolate within these we can pick a pair of opposing faces, draw a line perpendicular to them both through the point whose value we're approximating, use bilinear interpolation to approximate the values at the points of intersection of that line with each face and linearly interpolation between those points to approximate the point's value.
Finally, and I wouldn't be at all surprised if you've already guessed where we're going with this, in \(n\) dimensions the nodes sit at the corners of \(n\) dimensional rectangular cuboids which we can interpolate within by picking an opposing pair of \(n1\) dimensional hyperfaces, drawing a line perpendicular to them through the point, use \(n1\) dimensional multilinear interpolation to approximate the values at the points of intersection and finish the calculation with linear interpolation between them.
To understand why, first note that multilinear interpolation within an \(n\) dimensional rectangular cuboid is equivalent to a function
Secondly, note that the formula has \(2^n\) unknown coefficients which must be uniquely determined by the \(2^n\) nodes and values at the corners of the \(n\) dimensional cuboid and therefore the interpolation must be equivalent to exactly the same function no matter the order in which the hyperfaces are chosen.
Now you've probably noticed that the formula isn't actually linear but is instead an \(n^\mathrm{th}\) order polynomial. We can at least justify the name of the algorithm with the fact that it's conditionally linear in the sense that if we hold all but one of the \(x_i\) constant then it's linear in that one.
Here , where the subscripts \(i_j\) indicate that the node's \(i^\mathrm{th}\) element equals the \(j^\mathrm{th}\) tick on its axis which for a single cuboid will be either the lower or upper bound, indexed with zero and one respectively.
The values at the corner nodes of the opposing faces along the current axis are consequently stored in and and so this function simply linearly interpolates between the results of multilinear interpolation within those faces.
Program 1 uses it to interpolate within a multidimensional cuboid with randomly chosen lower and upper bounds having values at its corner nodes determined by a polynomial of the form that we proved is equivalent to multilinear interpolation and so should be exactly represented.
The differences between the interpolation and the polynomial are not inconsistent with those that we might expect from floating point rounding errors and we can therefore be confident that this recursive approach is correct.
For our
Listing 2 gives the
As it does so, it saves a copy of the grid axes in
It also keeps track of the indices of axes with more than two tick marks, together with a function for updating the guesses at their lower bounds of the cuboids that will contain the next point of interest, in the
The
Specifically, if there are three ticks on an axis then we can locate the lower bound by simply determining whether the element is less or greater than the second.
If there are more, but they are not too unevenly spaced so that
The
The indices of the tick marks upon each axis preceding the relevant elements of a point \(\mathrm{x}\) cannot be less than the first or greater than or equal to the last and so we can restrict the search for them to those within that range. Given that we can easily determine whether those elements lie before their current lower bounds or past their current upper bounds we can follow the decision as to whether we need to update them with a jump that restricts the new indices to just one or the other of those extremes as appropriate, as is done in listing 5.
Failing this, we can simply resort to a binary search through the sorted tick marks on each axis. Listing 6 uses our private by naming convention,
Here we are exploiting the fact that we know whether the tick that we're looking for is above or below the current lower bound to reduce the range of indices that we search through.
We'll be storing the values associated with the nodes in a nested array and so to initialise it with another nested array we simply need to recursively
When using a function of the vectors representing the node to calculate their associated values it's still easiest to use recursion construct the nested array of them, only this time we shall store the values of the tick marks as we iterate over them in an array at the index of the current axis, recursing to the next axis as we do so and finally using the fully populated array to construct an
Once we've initialised the values we need to check that they're consistent with the nodes and that they are all of the same type since we'll be implementing the interpolation for all of the arithmetic types in the
To support our own arithmetic types we'll need to use the
We shall slightly rearrange the formula for the interpolation from that implemented in listing 1 by passing down the accumulated products of the ratios of the differences between the elements and the tick marks preceding them and the differences between those enclosing them to cut down on the number of multiplication operations required for the calculation, as done by
The
The
The
Finally, you can find this implementation in MultiLinearInterpolate.js.
To give us a sense of what the results of multilinear interpolation looks like, program 5 uses
This illustration of the artefacts of multilinear interpolation, which are most clearly revealed along the curved valley floor, is as good a point as any to conclude.
[2] Cubic Line Division, www.thusspakeak.com, 2018.
[3] We're Not For Turning, www.thusspakeak.com, 2018.
[4] Cuboid Space Division, www.thusspakeak.com, 2018.
[5] I Still Haven't Found What I'm Looking For, www.thusspakeak.com, 2017.
Last time^{[4]}, in preparation for interpolating between multidimensional vector nodes, we implemented the
ak.grid
type to store ticks on a set of axes and map their intersections to ak.vector
objects to represent such nodes arranged at the corners of hyperdimensional rectangular cuboids.With this in place we're ready to take a look at one of the simplest multidimensional interpolation schemes; multilinear interpolation.
One Dimension, Two Dimension, Three Dimension, More
Figure 1: Rectangular Interpolation
To interpolate within one of those rectangles we can first choose a pair of opposing edges and draw a line perpendicular to them through the point of interest, then use linear interpolation along both of those edges to approximate the values at their points of intersection with that line and again between them to approximate the value at the point itself, as illustrated by figure 1 which marks the point of interest in red, the nodes in green and the points of intersection in black.
In three dimensions the nodes will sit at the corners of rectangular cuboids or, to put it more plainly, boxes. To interpolate within these we can pick a pair of opposing faces, draw a line perpendicular to them both through the point whose value we're approximating, use bilinear interpolation to approximate the values at the points of intersection of that line with each face and linearly interpolation between those points to approximate the point's value.
Finally, and I wouldn't be at all surprised if you've already guessed where we're going with this, in \(n\) dimensions the nodes sit at the corners of \(n\) dimensional rectangular cuboids which we can interpolate within by picking an opposing pair of \(n1\) dimensional hyperfaces, drawing a line perpendicular to them through the point, use \(n1\) dimensional multilinear interpolation to approximate the values at the points of intersection and finish the calculation with linear interpolation between them.
Any Which Way You Choose
One thing that might seem slightly dubious about this algorithm is that the order in which opposing hyperfaces are chosen isn't specified, but it turns out that it doesn't actually matter.To understand why, first note that multilinear interpolation within an \(n\) dimensional rectangular cuboid is equivalent to a function
\[
f_n\left(x_0,\dots,x_{n1}\right) = \sum_{i=0}^{2^n1} a_i \prod_{j=0}^{n1} x_j^{i_j}
\]
for some constants \(a_i\), where \(\sum\) is the summation sign, \(\prod\) is the product sign and \(i_j\) is the \(j^\mathrm{th}\) bit in the binary representation of \(i\) so that
\[
i = \sum_{j=0}^{n1} 2^{i_j}
\]
as proven in derivation 1.
Derivation 1: The Multilinear Interpolation Function
Firstly, since the structure of the formula is symmetric in \(x_i\), so that
Next, for some ordering of the axes of an \(n+1\) dimensional multilinear interpolation choose the opposing hyperfaces lying on the first \(n\) of them and draw a line through the point \(\left(x_0,\dots,x_{n1},x_n\right)\) to the points of intersection \(\left(x_0,\dots,x_{n1},x_{0,n}\right)\) and \(\left(x_0,\dots,x_{n1},x_{1,n}\right)\). If the formula is correct then their \(n\) dimensional multilinear interpolated values will be
\[
\sum_{i=0}^{2^n1} a_i \prod_{j=0}^{n1} x_j^{i_j} = \sum_{i=0}^{2^n1} a^\prime_i \prod_{j=0}^{n1} x_{p(j)}^{i_j}
\]
for any reordering \(p(j)\) of the indices \(j\), if it's correct for any ordering of the axes it must be correct for all of them.Next, for some ordering of the axes of an \(n+1\) dimensional multilinear interpolation choose the opposing hyperfaces lying on the first \(n\) of them and draw a line through the point \(\left(x_0,\dots,x_{n1},x_n\right)\) to the points of intersection \(\left(x_0,\dots,x_{n1},x_{0,n}\right)\) and \(\left(x_0,\dots,x_{n1},x_{1,n}\right)\). If the formula is correct then their \(n\) dimensional multilinear interpolated values will be
\[
\begin{align*}
f_{0,n}\left(x_0,\dots,x_{n1}\right) &= \sum_{i=0}^{2^n1} a_{0,i} \prod_{j=0}^{n1} x_j^{i_j}\\
f_{1,n}\left(x_0,\dots,x_{n1}\right) &= \sum_{i=0}^{2^n1} a_{1,i} \prod_{j=0}^{n1} x_j^{i_j}
\end{align*}
\]
for some constants \(a_{0,i}\) and \(a_{1,i}\). Linearly interpolating between these yields
\[
\begin{align*}
f_{n+1}(x_0,\dots,x_n) &= f_{0,n}\left(x_0,\dots,x_{n1}\right) + \frac{x_nx_{0,n}}{x_{1,n}x_{0,n}} \left(f_{1,n}\left(x_0,\dots,x_{n1}\right)f_{0,n}\left(x_0,\dots,x_{n1}\right)\right)\\
&= \left(\sum_{i=0}^{2^n1} a_{0,i} \prod_{j=0}^{n1} x_j^{i_j}\right) + \frac{x_nx_{0,n}}{x_{1,n}x_{0,n}} \left(\left(\sum_{i=0}^{2^n1} a_{1,i} \prod_{j=0}^{n1} x_j^{i_j}\right)  \left(\sum_{i=0}^{2^n1} a_{0,i} \prod_{j=0}^{n1} x_j^{i_j}\right)\right)\\
&= \left(\sum_{i=0}^{2^n1} a_{0,i} \prod_{j=0}^{n1} x_j^{i_j}\right) + \frac{x_nx_{0,n}}{x_{1,n}x_{0,n}} \sum_{i=0}^{2^n1} \left(a_{1,i}a_{0.i}\right) \prod_{j=0}^{n1} x_j^{i_j}\\
&= \left(\sum_{i=0}^{2^n1} \left(a_{0,i}\frac{a_{1,i}a_{0.i}}{x_{1,n}x_{0,n}}x_{0,n}\right) \prod_{j=0}^{n1} x_j^{i_j}\right) + \left(\sum_{i=0}^{2^n1}\frac{a_{1,i}a_{0.i} }{x_{1,n}x_{0,n}}x_n\prod_{j=0}^{n1} x_j^{i_j}\right)\\
&= \left(\sum_{i=0}^{2^n1} a^\prime_i x_n^0 \prod_{j=0}^{n1} x_j^{i_j}\right) + \left(\sum_{i=0}^{2^n1}a^\prime_{2^n+i} x_n^1\prod_{j=0}^{n1} x_j^{i_j}\right)
= \sum_{i=0}^{2^{n+1}1} a^\prime_i \prod_{j=0}^n x_j^{i_j}
\end{align*}
\]
so that the formula will be correct in \(n+1\) dimensions too. Finally, in one dimension we have
\[
f_1\left(x_0\right) = a_0 + a_1 x_0 = a_0 x_0^0 + a_1 x_0^1 = \sum_{i=0}^{2^11} a_i \prod_{j=0}^{11} x_j^{i_j}
\]
and so, by induction, it must be correct in all dimensions.
Secondly, note that the formula has \(2^n\) unknown coefficients which must be uniquely determined by the \(2^n\) nodes and values at the corners of the \(n\) dimensional cuboid and therefore the interpolation must be equivalent to exactly the same function no matter the order in which the hyperfaces are chosen.
Now you've probably noticed that the formula isn't actually linear but is instead an \(n^\mathrm{th}\) order polynomial. We can at least justify the name of the algorithm with the fact that it's conditionally linear in the sense that if we hold all but one of the \(x_i\) constant then it's linear in that one.
Implementing The Algorithm
Given the description of the multilinear interpolation algorithm, it should come as no surprise that the simplest way to implement it is by using recursion, as illustrated for a single hyperdimensional rectangular cuboid by listing 1.Listing 1: Recursive Multilinear Interpolation
function interpolate(i, x, lb, ub, values) { if(i===lb.length) return values; return (x[i]lb[i])/(ub[i]lb[i]) * interpolate(i+1, x, lb, ub, values[1]) + (ub[i]x[i])/(ub[i]lb[i]) * interpolate(i+1, x, lb, ub, values[0]); }
Here
i
is the index of the axis along which we're interpolating between faces, x
is an array representing the point whose value we're approximating, lb
and ub
are arrays of the lower and upper bounds of the cuboid upon each axis and values
is a nested array of the values at the corner nodes. For the latter we require that the value associated with a node \(\mathbf{x}_{i_0,i_1,\dots,i_{n1}}\) should be found at values[i_{0}][i_{1}]...[i_{n1}]
The values at the corner nodes of the opposing faces along the current axis are consequently stored in
values[0]
values[1]
Program 1 uses it to interpolate within a multidimensional cuboid with randomly chosen lower and upper bounds having values at its corner nodes determined by a polynomial of the form that we proved is equivalent to multilinear interpolation and so should be exactly represented.
Program 1: The Recursion Versus The Formula



The differences between the interpolation and the polynomial are not inconsistent with those that we might expect from floating point rounding errors and we can therefore be confident that this recursive approach is correct.
For our
ak
library's implementation we'll represent the nodes with an ak.grid
object, as stated in the introduction, and the values either as a nested array, as above, or use a function to calculate them from the ak.vector
objects representing the nodes that the former maps arrays of axis tick mark indices to.Listing 2 gives the
ak.multiLinearInterpolate
function which you can see first checks that the first argument is an ak.grid
object with at least two ticks per axis, all of which are finite, and then checks that the second argument is either an array or a function.Listing 2: ak.multiLinearInterpolate
ak.multiLinearInterpolate = function(grid, values) { var wide = []; var axes, dims, lb, dx, i, axis, n, t, interpolate, f; if(ak.type(grid)!==ak.GRID_T) { throw new Error('invalid grid in ak.multiLinearInterpolate'); } axes = grid.axes(); dims = axes.length; lb = new Array(dims); dx = new Array(dims); for(i=0;i<dims;++i) { axis = axes[i]; n = axis.length; if(n<2) { throw new Error('axis too short in ak.multiLinearInterpolate'); } if(!isFinite(axis[0])  !isFinite(axis[n1])) { throw new Error('nonfinite axis in ak.multiLinearInterpolate'); } if(n>2) wide.push({i:i, lb:updateLB(axes[i])}); lb[i] = 0; dx[i] = axis[1]axis[0]; } switch(ak.nativeType(values)) { case ak.ARRAY_T: values = copyValues(values); break; case ak.FUNCTION_T: values = mapValues(0, dims, new Array(dims), axes, values); break; default: throw new Error('invalid values in ak.multiLinearInterpolate'); } t = valueType(values); validateValues(0, dims, t, grid, values); interpolate = interpolateType(t, dims, values); f = function(x) {return interpolate(transformX(x, wide, axes, lb, dx), lb);}; f.grid = function() {return grid;}; f.values = function() {return copyValues(values);}; f.nodes = function() {return makeNodes(dims, grid, values);}; return Object.freeze(f); };
As it does so, it saves a copy of the grid axes in
axes
, the dimension of the interpolation in dims
and initialises the lb
and dx
arrays with a first guess of zero at the indices of the tick marks that will immediately precede the first point whose value we shall approximate and the distances between them and those following them.It also keeps track of the indices of axes with more than two tick marks, together with a function for updating the guesses at their lower bounds of the cuboids that will contain the next point of interest, in the
wide
array, initialises the values
with the copyValues
helper function if they are passed as a nested array and with mapValues
if a function is given to calculate them, validates the values and defines an interpolate
function based upon their type before finally creating a function to perform the interpolation and assign it property access methods to return the grid defining the nodes, the nested array of values and the set of nodes themselves.The
updateLB
helper function uses exactly the same tricks as our one dimensional, or univariate, interpolation algorithms used to efficiently locate the tick marks immediately preceding each element of a vector representing a point upon its associated axis, as illustrated by listing 3.Listing 3: The updateLB Helper Function
function updateLB(axis) { var x0, dx; if(axis.length===3) { return function(x) {return x<axis[1]?0:1;}; } x0 = axis[0]; dx = axis[axis.length1]x0; if(canJump(axis, x0, dx)) { return function(x, lb) {return jumpLB(x, axis, x0, dx, lb);}; } return function(x, lb) {return findLB(x, axis, lb);}; }
Specifically, if there are three ticks on an axis then we can locate the lower bound by simply determining whether the element is less or greater than the second.
If there are more, but they are not too unevenly spaced so that
\[
x_{i1} \leqslant \bigg\lfloor\frac{xx_0}{x_{n1}x_0} \times n\bigg\rfloor \leqslant x_{i+1}
\]
for
\[
x_i \leqslant x \leqslant x_{i+1}
\]
where the oddly shaped brackets mean the largest integer that is less than or equal to the term between them, then we can use
\[
\bigg\lfloor\frac{xx_0}{x_{n1}x_0} \times n\bigg\rfloor
\]
as a first guess at the index of the tick immediately preceding \(x\) and be sure that in the worst case we'll only be off by one in either direction.The
canJump
function given in listing 4 determines whether or not this is the case by checking it for the tick marks themselves, since this will mean that it must be for any \(x\) between them, ignoring the first and last because they must be identified as the lower and upper bounds of the first and last pair respectively.Listing 4: The canJump Helper Function
function canJump(axis, x0, dx) { var n = axis.length; var i, j; for(i=1;i<n1;++i) { j = ak.floor(n*((axis[i]x0)/dx)); if(j<i1  j>i) return false; } return true; }
The indices of the tick marks upon each axis preceding the relevant elements of a point \(\mathrm{x}\) cannot be less than the first or greater than or equal to the last and so we can restrict the search for them to those within that range. Given that we can easily determine whether those elements lie before their current lower bounds or past their current upper bounds we can follow the decision as to whether we need to update them with a jump that restricts the new indices to just one or the other of those extremes as appropriate, as is done in listing 5.
Listing 5: Jumping To The Lower Bound
function jumpLBFwd(x, axis, x0, dx) { var n = axis.length; var i = Math.min(ak.floor(n*((xx0)/dx)), n2); if(x<axis[i]) {if(i>0) i;} else if(x>axis[i+1]) {if(i<n2) ++i;} return i; } function jumpLBRev(x, axis, x0, dx) { var n = axis.length; var i = Math.max(ak.floor(n*((xx0)/dx)), 0); if(x<axis[i]) {if(i>0) i;} else if(x>axis[i+1]) {if(i<n2) ++i;} return i; } function jumpLB(x, axis, x0, dx, lb) { if(x>axis[lb+1]) { if(lb<axis.length2) lb = jumpLBFwd(x, axis, x0, dx); } else if(x<axis[lb]) { if(lb>0) lb = jumpLBRev(x, axis, x0, dx); } return lb; }
Failing this, we can simply resort to a binary search through the sorted tick marks on each axis. Listing 6 uses our private by naming convention,
ak._unsafeLowerBound
function^{[5]} to do so since we know that we are already satisfying the conditions that the public ak.lowerBound
would have checked for.Listing 6: Searching For The Lower Bound
function compare(l, r) {return lr;} function findLBFwd(x, axis, lb) { if(++lb<axis.length2 && x>axis[lb+1]) { lb = ak._unsafeLowerBound(axis, x, compare, lb+2, axis.length1)1; } return lb; } function findLBRev(x, axis, lb) { if(lb>0 && x<axis[lb]) { lb = ak._unsafeLowerBound(axis, x, compare, 1, lb)1; } return lb; } function findLB(x, axis, lb) { if(x>axis[lb+1]) { if(lb<axis.length2) lb = findLBFwd(x, axis, lb); } else if(x<axis[lb]) { if(lb>0) lb = findLBRev(x, axis, lb); } return lb; }
Here we are exploiting the fact that we know whether the tick that we're looking for is above or below the current lower bound to reduce the range of indices that we search through.
We'll be storing the values associated with the nodes in a nested array and so to initialise it with another nested array we simply need to recursively
slice
it, its elements, its elements' elements and so on until we get to a nonarray element that we can directly copy, as illustrated by listing 7.Listing 7: Recursively Copying The Values
function copyValues(values) { var n; if(ak.nativeType(values)===ak.ARRAY_T) { n = values.length; values = values.slice(0); while(n>0) values[n] = copyValues(values[n]); } return values; }
When using a function of the vectors representing the node to calculate their associated values it's still easiest to use recursion construct the nested array of them, only this time we shall store the values of the tick marks as we iterate over them in an array at the index of the current axis, recursing to the next axis as we do so and finally using the fully populated array to construct an
ak.vector
object to pass to the function once we've run out of axes, as shown by listing 8.Listing 8: Recursively Calculating The Values
function mapValues(i, dims, x, axes, f) { var axis, n, values; if(i<dims) { axis = axes[i]; n = axis.length; values = new Array(n); while(n>0) { x[i] = axis[n]; values[n] = mapValues(i+1, dims, x, axes, f); } } else { values = f(ak.vector(x)); } return values; }
Once we've initialised the values we need to check that they're consistent with the nodes and that they are all of the same type since we'll be implementing the interpolation for all of the arithmetic types in the
ak
library that support addition with themselves and multiplication by numbers, just as we did for ak.linearInterpolate
. The valueType
function given in listing 9 recursively extracts the first element of each of the nested arrays and finally determines the type of the values it contains with ak.type
. The validateValues
function recursively checks that each nested array has the same length as the number of tick marks upon their associated axis and that they contain values of type t
.Listing 9: Validating The Values
function valueType(values) { return (ak.nativeType(values)===ak.ARRAY_T) ? valueType(values[0]) : ak.type(values); } function validateValues(i, dims, t, grid, values) { var n; if(i<dims) { n = grid.length(i); if(ak.nativeType(values)!==ak.ARRAY_T  values.length!==n) { throw new Error('values size mismatch in ak.multiLinearInterpolate'); } while(n>0) validateValues(i+1, dims, t, grid, values[n]); } else { if(ak.type(values)!==t) { throw new Error('values type mismatch in ak.multiLinearInterpolate'); } } }
To support our own arithmetic types we'll need to use the
ak.add
and ak.mul
overloads to calculate the interpolated values but since this will be inefficient if the values are numbers it makes sense to treat them as a special case using JavaScript's built in operators. The interpolateType
function given in listing 10 switches between using the interpolateNumber
function for numbers and the interpolateGeneral
function for types that have appropriate overloads of ak.add
and ak.mul
.Listing 10: Choosing An Implementation
function interpolateType(t, dims, values) { var add, mul; if(t===ak.NUMBER_T) { return function(x, lb) { return interpolateNumber(1, 0, dims, values, x, lb); }; } try{add=ak.add[t][t]; mul=ak.mul[ak.NUMBER_T][t];} catch(e){} if(ak.nativeType(add)!==ak.FUNCTION_T  ak.nativeType(mul)!==ak.FUNCTION_T) { throw new Error('nonarithmetic value type in ak.multiLinearInterpolate'); } return function(x, lb) { return interpolateGeneral(add, mul, 1, 0, dims, values, x, lb); }; }
We shall slightly rearrange the formula for the interpolation from that implemented in listing 1 by passing down the accumulated products of the ratios of the differences between the elements and the tick marks preceding them and the differences between those enclosing them to cut down on the number of multiplication operations required for the calculation, as done by
interpolateNumber
in listing 11.Listing 11: Interpolating Numbers
function interpolateNumber(vol, i, dims, values, x, lb) { var xi, j; if(i===dims) return vol*values; xi = x[i]; j = lb[i]; return interpolateNumber(vol*xi[0], i+1, dims, values[j+1], x, lb) + interpolateNumber(vol*xi[1], i+1, dims, values[j], x, lb); }
The
interpolateGeneral
function given in listing 12 performs the same calculation using our ak.add
and ak.mul
operators for general arithmetic types.Listing 12: Interpolating General Arithmetic Types
function interpolateGeneral(add, mul, vol, i, dims, values, x, lb) { var xi, j; if(i===dims) return mul(vol, values); xi = x[i]; j = lb[i]; return add(interpolateGeneral(add, mul, vol*xi[0], i+1, dims, values[j+1], x, lb), interpolateGeneral(add, mul, vol*xi[1], i+1, dims, values[j], x, lb)); }
The
transformX
function given in listing 13 first checks that the argument of the interpolation is an ak.vector
object with the correct number of dimensions, then updates the lower bounds on each axis with more than two tick marks using the lb
function that we chose as being the most efficient and, if it has changed, updates it and the width of the interval between it and the next. It then replaces the elements of the array converted from that argument with the pair of ratios of the differences between them and the lower and upper tick marks with the differences between those lower and upper tick marks.Listing 13: The transformX Helper Function
function transformX(x, wide, axes, lb, dx) { var n, i, axis, lbi, dxi, xi; if(ak.type(x)!==ak.VECTOR_T) { throw new Error('invalid argument in ak.multiLinearInterpolate'); } if(x.dims()!==axes.length) { throw new Error('argument size mismatch in ak.multiLinearInterpolate'); } x = x.toArray(); n = wide.length; while(n>0) { i = wide[n].i; axis = axes[i]; lbi = wide[n].lb(x[i], lb[i]); if(lbi!==lb[i]) { lb[i] = lbi; dx[i] = axis[lbi+1]axis[lbi]; } } i = x.length; while(i>0) { axis = axes[i]; lbi = lb[i]; dxi = dx[i]; xi = x[i]; x[i] = [(xiaxis[lbi])/dxi, (axis[lbi+1]xi)/dxi]; } return x; }
The
makeNodes
function given in listing 14 converts the grid and nested array of values into a set of node/value pairs using the former's facility for mapping from arrays of indices to ak.vector
objects and the ak.nextState
function to iterate over them.Listing 14: Converting The Grid And Values Into Nodes And Values
function makeNodes(dims, grid, values) { var ub = grid.lengths(); var pos = new Array(dims); var n = 1; var i, nodes, y; for(i=0;i<dims;++i) {pos[i]=0; n*=ub[i];} nodes = new Array(n); n = 0; do { y = values; for(i=0;i<dims;++i) y = y[pos[i]]; nodes[n++] = {x:grid.map(pos), y:y}; } while(ak.nextState(pos, ub)); return nodes; }
Finally, you can find this implementation in MultiLinearInterpolate.js.
Using The Interpolator
Program 2 demonstratesak.multiLinearInterpolate
by using it to interpolate between nodes with associated values calculated by the formula representing the results of multilinear interpolation, just as we did in program 1.
Program 2: The Interpolator Versus The Formula



To give us a sense of what the results of multilinear interpolation looks like, program 5 uses
ak.multiLinearInterpolate
to approximate Rosenbrock's function
\[
f\left(x_0, x_1\right) = 100 \times \left(x_1x_0^2\right)^2 + \left(1x_0\right)^2
\]
plotting the function itself at the top and the interpolation beneath.
Program 3: Interpolating Rosenbrock's Function



This illustration of the artefacts of multilinear interpolation, which are most clearly revealed along the curved valley floor, is as good a point as any to conclude.
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.
[4] Cuboid Space Division, www.thusspakeak.com, 2018.
[5] I Still Haven't Found What I'm Looking For, www.thusspakeak.com, 2017.
Leave a comment