⎕io=0 assumed throughout; works in 1-origin with the obvious modifications.
Introduction
On Wednesday, a question arrived via Dyalog Support from an intern in Africa: If M is the matrix on the left, use linear interpolation to compute the result on the right.
1 20 1 20
4 80 2 40
6 82 3 60
4 80
5 81
6 82
Linear Interpolation
Two points (x0,y0) and (x1,y1) specify a line; for any x there is a unique y on that line (assuming x0≠x1). The equation for the line derives as follows, starting from its slope m:
m = (y1-y0) ÷ (x1-x0)
(y-y0) = m × (x-x0)
y = y0 + m × (x-x0)
Therefore, if ⍺ is a 2-by-2 matrix of the two points and ⍵ are the x-values to be interpolated, then:
g ← {(⊃⌽⍺)+(⍵-⊃⍺)÷÷/-⌿⍺}
⊢ M←1 4 6,⍪20 80 82
1 20
4 80
6 82
M[0 1;] g 2 3
40 60
M[1 2;] g 5
81
A New Twist, A New Solution
The problem as posed implicitly required that:
- The x-values are the positive integers bounded by
⊃⊖M. - Appropriate rows of the matrix are selected for a given x-value.
- The missing x-values and their interpolations are “slotted back” into the argument matrix.
These requirements are best met by ⍸, interval index, a relatively new primitive function introduced in Dyalog APL version 16.0. The left argument ⍺ must be sorted and partitions the universe into disjoint contiguous intervals; ⍺⍸⍵ finds the index of the interval which contains an item of ⍵. The result is ⎕io dependent.
For the given matrix M, the partition (of the real numbers in this case) is depicted below. As in conventional mathematical notation, [ denotes that the interval includes the left end-point and ) denotes that the interval excludes the right end-point.
1 4 6
─────────)[───────)[─────)[──────────
¯1 0 1 2
v←¯5 0 1 2.5 6 3 4 5 9 8 7
1 4 6 ⍸ v
¯1 ¯1 0 0 2 0 1 1 2 2 2
v ,[¯0.5] 1 4 6 ⍸ v
¯5 0 1 2.5 6 3 4 5 9 8 7
¯1 ¯1 0 0 2 0 1 1 2 2 2
With ⍸ in hand, the problem can be solved as follows:
interpol←{
(x y)←↓⍉⍵
m←m,⊃⌽m←(2-/y)÷(2-/x)
j←0⌈x⍸i←1+⍳⊃⌽x
i,⍪y[j]+m[j]×i-x[j]
}
interpol M
1 20
2 40
3 60
4 80
5 81
6 82
The problem of x-values less than the first end-point is finessed by applying 0⌈ to the interval indices, and that of x-values greater than or equal to the last end-point is finessed by repeating the last slope m←m,⊃⌽m.
It is possible to do the interpolation only on the missing indices (2 3 5 in this case) and insert them into the argument matrix. It seems neater to simply interpolate everything, and in so doing provide a check that the interpolated values equal the values given in the argument.
An Alternative Interpolation
Interpolating according to two selected rows of a matrix of points treats the function as piecewise linear, with sharp inflection points where the lines join (different slopes between adjacent lines). A “holistic” alternative approach is possible: the matrix can be interpreted as specifying a single line and the interpolation is according to this single line. The ⌹ primitive function computes the coefficients of the line which best fits the points:
⎕rl←7*5 ⍝ for reproducible random numbers
⊢ M←t,⍪(?7⍴5)+¯17+3×t←?7⍴100
35 89
98 278
19 44
4 ¯5
62 170
49 133
25 59
M[;1] ⌹ 1,M[;,0] ⍝ y-intercept and slope
¯15.3164 2.99731
interpola ← {(1,⍤0⊢⍵)+.×⍺[;1]⌹1,⍺[;,0]}
M[;1] ,[¯0.5] M interpola M[;0]
89 278 44 ¯5 170 133 59
89.5895 278.42 41.6325 ¯3.32713 170.517 131.552 59.6164
M interpola 33 35 37 39.7
83.5949 89.5895 95.5841 103.677
Finally
Our best wishes to the intern. Welcome to APL!
