New Life for Old Code

Robert B. Nachbar, Wolfram Solutions

ABSTRACT

A lot of good code was written in FORTRAN in the days of punch cards. Unfortunately, a good deal of it was lost along with the boxes of cards. I recently had need of Akima’s spline method that we used a lot in those days, and reimplementing the algorithm from the original FORTRAN (still online at acm.org) in Wolfram Language was less than ideal. Since then I’ve had the opportunity to start fresh from the original math. This talk will describe the new implementation of both the spline and interpolation methods as CompiledFunction and InterpolatingFunction objects, and their deployment in the Wolfram Function Repository. Additionally, Akima had published a 2D interpolation algorithm, and I will describe how his method can be extended to n-dimensions and implemented efficiently with ListConvolve.

Interpolation

From the Meriam-Webster Dictionary

interpolation noun
in·​ter·​po·​la·​tion | \ in-ˌtər-pə-ˈlā-shən

a : an act of interpolating something or the state of being interpolated : the introduction or insertion of something spurious or foreign

… versions disfigured by the frequent and substantial interpolation of freely invented matter …

— Bernard Knox

b : something that is introduced or inserted : an insertion or addition

How can such interpolations be identified, and what literary processes led to these additions?

— Jan Nattier

: the process of calculating an approximate value based on values that are already known

The age of tephras that were not directly dated was estimated by linear interpolation from the calibrated date.

— Richard J. Payne and Jeffrey J. Blackford

From Wikipedia

In the mathematical field of numerical analysis, interpolation is a type of estimation, a method of constructing new data points within the range of a discrete set of known data points.

In Wolfram Language

\[KernelIcon];?Interpolation

\[KernelIcon];?FindFit

Interpolation is a piecewise continuous method

passes through each point of the input

typically a polynomial of degree n is used and fit to n+1 points about each interval

“smoothness” controlled degree of polynomial

Fitting presumes a specific functional form that spans the domain

usually passes close to each point of the input

\[KernelIcon];example[1]

Akima’s Method for a Line

Interpolation

The smoothness of an interpolation is controlled by the order or degree of the polynomial

often fails when the set of points have “sharp corners”

increasing the order helps smooth the curve, but then can “over shoot” extrema

\[KernelIcon];example[2]

Akima’s interpolation

Akima’s method smooths the interpolation by providing slopes along with the  data values

A cubic polynomial is fit to each interval

a+b x+c (x^2)+d (x^3) ⟶ 4 unkonwns

({{(Subscript[x, 1]), (Subscript[y, 1])}, {(Subscript[x, 2]), (Subscript[y, 2])}, {(Subscript[x, 1]), (Subscript[s, 1])}, {(Subscript[x, 2]), (Subscript[s, 2])}})

the interpolation is continuous and smooth

\[KernelIcon];example[3]

The FORTRAN subroutines were published in Communications ACM

We punched them up (yes, punch cards!) and used them for many chemistry and physics labs

The Antecedent (FORTRAN)

Akima, H. “A new method of interpolation and smooth curve fitting based on local procedures” J. Assoc. Comput. Machinery, 17 (4), 589-602 (1970).

Akima, H. “Algorithm 433. Interpolation and smooth curve fitting based on local procedures [E2]” Commun. ACM, 15 (10), 914-918 (1972).

\[KernelIcon];theCode[]

One of the drawbacks of the original algorithm was that the user had to provide the points at which the interpolation was desired

Sometimes this led to plots that did not appear smooth

Long story short

The back story

Going back to the original publication and learning the math:

“With five data points 1, 2, 3, 4, and 5, as shown in Figure 3, the author [2] previously proposed a conditon that the slope of line segment (OverBar[CD]), tangent to the curve at point 3, is determined by

|(OverBar[(2C)])/(OverBar[CA])|==|(OverBar[4D])/(OverBar[DB])|

where the point of intersection of the two straight lines extended from line segments (OverBar[12]) and (OverBar[34]) is denoted by A, …”

Which can be implemented very susinctly in WL:

intervalSlopes = Divide @@@ Subtract @@@ Partition[Reverse /@ data, 2, 1]; 
		
pointSlopes = With[{d = Abs[#4 - #3] + Abs[#2 - #1]}, 
	If[d != 0, (Abs[#4 - #3] #2 + Abs[#2 - #1] #3)/d, (#2 + #3)/2]] & @@@ 
	Partition[intervalSlopes, 4, 1]; 

Interpolation[MapThread[{Most[#1], Last[#1], #2} &, {data[[3 ;; -3]], pointSlopes}], 
	PeriodicInterpolation -> periodic]

Implementation details

The data need to be extrapolated by 2 points on each end of the data range

Akima provided a useful heuristic

The points are used only for computation of slopes at the first two and last two points

The WL function Interpolation takes not only ((Subscript[x, i]),(Subscript[y, i])) data, but also derivatives ((Subscript[x, i]),(Subscript[y, i]^′)) data

The resulting InterpolatingFunction computes the exact same results as the FORTRAN code

It’s efficient, differentiable, continuous, and can be used like any other numeric function

The Wolfram Function Repository (WFR)

The WFR is a convenient way of publishing your code

It’s easy to make your efforts look polished and professional

As a ResourceFunction, one does not need to load a package

\[KernelIcon];ro=ResourceObject["AkimaInterpolation"]

\[KernelIcon];ro["Description"]

Here is the reference page:

\[KernelIcon];ro["DocumentationNotebook"];

Using the function

\[KernelIcon];ResourceFunction["AkimaInterpolation"]

\[KernelIcon];data={{0,0},{1,0},{2,0},{3,0},{4,0},{5,1},{6,10},{7,80},{8,100},{9,150}};
afun=({{"[◼]", "AkimaInterpolation"}})[data];
ifun=Interpolation[data];
sfun=Interpolation[data,Method->"Spline"];

\[KernelIcon];Plot[{ifun[x],sfun[x],afun[x]},{x,0,9},Epilog->{Red,PointSize[0.02],Point[data]},PlotLegends->{"Hermite","Spline","Akima"},AspectRatio->1]

Zoom in to the problem area

\[KernelIcon];Plot[{ifun[x],sfun[x],afun[x]},{x,3.75,6.25},Epilog->{Red,PointSize[0.02],Point[data]},PlotLegends->{"Hermite","Spline","Akima"},AspectRatio->1]

Periodic data

\[KernelIcon];pfun=({{"[◼]", "AkimaInterpolation"}})[{1,5,7,2,3,1},PeriodicInterpolation->True]

\[KernelIcon];Plot[pfun[x],{x,-12,18}]

Akima’s Spline

Akima also published a spline interpolation, one that allows multiple ordinate values at the same abscissa value

\[KernelIcon];ResourceFunction["AkimaSpline"]

\[KernelIcon];data={{-30,70},{-30,40},{-30,20},{-18,4},{0,0},{18,4},{30,20},{30,40},{30,50},{30,70}};

\[KernelIcon];sfun=({{"[◼]", "AkimaSpline"}})[data]

In this case, a CompiledFunction is returned that computes a parametic result

\[KernelIcon];ParametricPlot[sfun[z],{z,0,1},Axes->False,Frame->True,AspectRatio->1,Epilog->{Red,PointSize[0.02],Point[data]}]

A more interesting example

\[KernelIcon];data={{1,1},{0,-0.5},{-1,1},{0.5,0},{-1,-1},{0,0.5},{1,-1},{-0.5,0},{1,1}};
ifun=({{"[◼]", "AkimaSpline"}})[data,PeriodicInterpolation->True];
ParametricPlot[ifun[z],{z,0,1},Axes->False,Frame->True,AspectRatio->1,PlotRange->All,Epilog->{Red,PointSize[0.02],Point[data]}]

Motivation for this talk

We wanted to do something outrageous for the Wolfram Summer School

First, get the data that define a particular graph

\[KernelIcon];spikeyData=;

\[KernelIcon];points=FirstCase[spikeyData,GraphicsComplex[p_,_]:>p,{},∞];
edges=UndirectedEdge@@@Join@@(Partition[#,2,1]&)/@Cases[spikeyData,Line[p_]:>p,∞];
gr=Graph[edges,ImageSize->Small]

Then, find all the cycles with 6 vertices

\[KernelIcon];cycles=FindCycle[gr,{6},All];
vertexLists=Function[cycle,With[{v=First/@cycle},Append[v,First@v]]]/@cycles;
lists=With[{keyVertices={3,9,15,21,27}},Cases[Cases[vertexLists,{___,#,___}],l_/;FreeQ[l,Alternatives@@DeleteCases[keyVertices,#]]]&/@keyVertices];

Now, make the parametric plots

\[KernelIcon];curves=({{"[◼]", "AkimaSpline"}})[N@points[[#]],PeriodicInterpolation->True]&/@Join@@lists;ParametricPlot[Through[curves[u]],{u,0,1},Axes->False,Frame->False]

Lastly, leave one’s mark on the world

Akima’s Method for a Surface

Akima also published an extension of his method to 2-dimensional surfaces:

Akima, H. “A method of bivariate interpolation and smooth surface fitting based on local procedures” Commun. ACM, 17 (1), 18-20 (1974).

Akima, H. “Algorithm 474. Bivariate interpolation and smooth surface fitting based on local procedures [E2]” Commun. ACM, 17 (1), 26-31 (1974).

This WL code is not in WFR yet, so load the local package

\[KernelIcon];<<Akima`

Example

Import some data

\[KernelIcon];{dataRange2D,data2D}=With[{d=Transpose@N@Import[FindFile["Akima74bivariate--data.csv"]]},{{Rest@First@d,First/@Rest@d},Rest/@Rest@d}];
Dimensions@data2D
Length/@dataRange2D

\[KernelIcon];ListPlot3D[data2D,DataRange->MinMax/@dataRange2D,BoxRatios->Automatic]

\[KernelIcon];?AkimaInterpolation2D

\[KernelIcon];surface=AkimaInterpolation2D[data2D,DataRange->dataRange2D]

\[KernelIcon];Row[{Plot3D[surface[x,y],{x,0,50},{y,0,40},BoxRatios->Automatic],ContourPlot[surface[x,y],{x,0,50},{y,0,40},AspectRatio->Automatic,Contours->{-0.1,0.1,1,5,10,20,30,40,50,60}]},Spacer[36]]

Implementation

Again, the heart of the algorithm is compact

(* compute interval gradients *)
c = Transpose @ MapThread[Divide, {Subtract @@@ Partition[Transpose @ z, 2, 1], Subtract @@@ Partition[x, 2, 1]}];
d = MapThread[Divide, {Subtract @@@ Partition[z, 2, 1], Subtract @@@ Partition[y, 2, 1]}]; 
e = MapThread[Divide, {Subtract @@@ Partition[c, 2, 1], Subtract @@@ Partition[y, 2, 1]}];
		
wx = Subtract @@@ Partition[Transpose @ c, 2, 1] // Transpose // Abs; 
wy = Subtract @@@ Partition[d, 2, 1] // Abs; 
		
(* compute derivatives *)
zx = Function[cj, With[{w2 = Abs[#4 - #3], w3 = Abs[#2 - #1]}, 
	If[w2 + w3 != 0, (w2 #2 + w3 #3)/(w2 + w3), (#2 + #3)/2]
	] & @@@ Partition[cj, 4, 1]] /@ c[[3 ;; -3]];
zy = Function[di, With[{w2 = Abs[#4 - #3], w3 = Abs[#2 - #1]}, 
	If[w2 + w3 != 0, (w2 #2 + w3 #3)/(w2 + w3), (#2 + #3)/2]
	] & @@@ Partition[di, 4, 1]] /@ Transpose[d][[3 ;; -3]] // Transpose;
zxy = Table[
	(wx[[i, j]] (wy[[i, j]] e[[i - 1, j - 1]] + wy[[i - 2, j]] e[[i, j - 1]]) + 
		wx[[i, j - 2]] (wy[[i, j]] e[[i - 1, j]] + wy[[i - 2, j]] e[[i, j]])) / 
	(If[wx[[i, j]] + wx[[i, j - 2]] == 0., 1., wx[[i, j]] + wx[[i, j - 2]]] * 
		If[wy[[i, j]] + wy[[i - 2, j]] == 0., 1., wy[[i, j]] + wy[[i - 2, j]]]), 
	{i, 3, Length[y] - 2}, 
	{j, 3, Length[x] - 2}
	];
		
(* interpolate *)
Interpolation[Join @@ Table[{{xIn[[j]], yIn[[i]]}, zIn[[i, j]], {zx[[i, j]], zy[[i, j]]}, 
	{{Automatic, zxy[[i, j]]}, {zxy[[i, j]], Automatic}}}, {i, 1, Length[yIn]}, {j, 1, Length[xIn]}]]

Extension to n Dimenions

As we saw for the 1-D and 2-D cases, the method for the interval slopes, or divided differences, is the same, as are the first derivatives ((dy)/(dx)) and ((dy)/(dx))

Second order divided differences are computed from first order divided differences, and second derivatives (((d^2)y)/(dxdy)) are computed with a bilinear form

These observations suggest that this method should be ammenable to extension to n-dimensions

This is still a work in progress!

Fitting higher dimension cubic polynomials

Don’t have a proof, but have demonstrated symbolic solutions for the coefficients up through n==3 using just function values, first derivatives, and “polylinear” derivatives at the points of the hypercube

Computing interval slopes

To compute the interval gradients for the n==2 case, we used a conjugation operation to allow reuse of the code from the x-dimension for the y-dimension.

zx = Function[cj, With[{w2 = Abs[#4 - #3], w3 = Abs[#2 - #1]}, 
	If[w2 + w3 != 0, (w2 #2 + w3 #3)/(w2 + w3), (#2 + #3)/2]
	] & @@@ Partition[cj, 4, 1]] /@ c[[3 ;; -3]];
zy = Function[di, With[{w2 = Abs[#4 - #3], w3 = Abs[#2 - #1]}, 
	If[w2 + w3 != 0, (w2 #2 + w3 #3)/(w2 + w3), (#2 + #3)/2]
	] & @@@ Partition[di, 4, 1]] /@ Transpose[d][[3 ;; -3]] // Transpose;

For higher values of n, this gets unwieldy

Instead, we can use ListConvolve

Think edge detection in images

We just need a difference kernel for each dimension

\[KernelIcon];generateKernels[kernel_,n_] := With[{basis =Nest[List, kernel, n - 1], perms =Table[SubsetMap[Reverse, Range[n], ij], {ij, Reverse@Partition[Range[n], 2, 1]}]},
  FoldList[Transpose, basis, perms]
  ]

\[KernelIcon];generateKernels[{1,-1},#]&/@Range[4]//

Generate some 4-dimensional data

\[KernelIcon];data=RandomInteger[{-5,10},{4,3,5,7}];
Dimensions@%

Generate the values along the grid edges

\[KernelIcon];dataRange=Accumulate@RandomInteger[{1,5},#]&/@Dimensions[data]

\[KernelIcon];n=Length@Dimensions[data]

Generate the kernels, one for each dimension

\[KernelIcon];kernels=generateKernels[{1,-1},n]

Generate the denominators for the interval derivatives

\[KernelIcon];deltas=ListConvolve[{1,-1},#]&/@dataRange

Generate the levels at whcih to apply the kernels

\[KernelIcon];levels=List/@Range[0,n-1]

Compute the (1^st) degree divided differences (interval derivatives)

\[KernelIcon];df=MapThread[Function[{ker,del,lev},Map[MapThread[Divide,{#,del}]&,ListConvolve[ker,data,{-1,1},{},Times,Plus,n],lev]],{generateKernels[{1,-1},n],Reverse@deltas,Reverse@levels}];
Dimensions/@%

(2^nd) and higher degree and higher divided differences are computed from lower degree divided differences

We can simulate the process to see how well it works

\[KernelIcon];dFlabels[base_String, kernels : {__String}] := Module[{dy = StringJoin[base, #] & /@ kernels},
  Join[dy, MapThread[dFlabels, {Most@dy, Rest@NestList[Rest, kernels, Length@kernels - 1]}]]
  ]

\[KernelIcon];ddLabels=dFlabels["",{"w","x","y","z"}]

However, they are in a somewhat cumbersome order

Simply flattening the nested lists is not useful, but we can extract level by level

\[KernelIcon];Function[r,Join@@Echo[(Cases[r,_String,{#1}]&)/@Range[Depth[r]-1]]][ddLabels]

Now, for the actual data

\[KernelIcon];dividedDifferences[f_,kernels_,deltas_,levels_]:=
Module[{df,n},
n=Length@Dimensions@f;
df=MapThread[Function[{ker,del,lev},Map[MapThread[Divide,{#,del}]&,ListConvolve[ker,f,{-1,1},{},Times,Plus,n],lev]],{kernels,deltas,levels}];
Join[df,MapThread[dividedDifferences,{Most@df,Rest@NestList[Rest,kernels,Length@kernels-1],Rest@NestList[Rest,deltas,Length@deltas-1],Rest@NestList[Rest,levels,Length@levels-1]}]]
]

\[KernelIcon];extractByDegree[tree_]:=With[{d=Depth@First@tree},Join@@(Cases[tree,_?(Depth[#]==d&),{#}]&)/@Range@(Depth@tree-1)]

\[KernelIcon];ddData=dividedDifferences[data,kernels,Reverse@deltas,Reverse@levels]//extractByDegree;
Dimensions/@%

computing derivatives

The steps required to complete the process have not been fully worked out

However, it is instructive to show how ListConvolve can be used to select the non-adjacent elements from the weight matrices needed to compute the final derivatives

\[KernelIcon];(data=Array[(Subscript[w, ##])&,{5,5}])//MatrixForm

\[KernelIcon];example[4]

\[KernelIcon];n=Length@Dimensions@data

\[KernelIcon];kernels=generateKernels[{True,False,True},2]

\[KernelIcon];Grid/@Map[Apply[Join,ListConvolve[#1,data,{-1,1},{},Pick[#2,#1]&,List,n],{n}]&,kernels]

Thank you!

Initialization

SetDirectory@NotebookDirectory[];

PrependTo[$Path,"/Users/rnachbar/Documents/Software/WolframResearch/src/Akima/"];

example[1]:=DynamicModule[{data={{0,0},{1,0.8415},{2,0.9093},{3,0.1411},{4,-0.7568},{5,-0.9589},{6,-0.2794}},hfun1,sfun,hfun3},
hfun1=Interpolation[data,InterpolationOrder->1];
hfun3=Interpolation[data,InterpolationOrder->3];
sfun=Interpolation[data,Method->"Spline"];

TabView[{"  data  "->ListPlot[data,PlotStyle->Directive[Red],PlotMarkers->{Automatic,Large},GridLines->{None,{0}},ImageSize->Large],Column[{"linear","interpolation"},Center]->Show[ListPlot[data,PlotStyle->Directive[Red],PlotMarkers->{Automatic,Large},GridLines->{None,{0}},ImageSize->Large],Plot[hfun1[x],{x,0,6},PlotStyle->Directive[Thick],GridLines->{None,{0}},ImageSize->Large]],Column[{"cubic","interpolation"},Center]->Show[ListPlot[data,PlotStyle->Directive[Red],PlotMarkers->{Automatic,Large},GridLines->{None,{0}},ImageSize->Large],Plot[hfun3[x],{x,0,6},PlotStyle->Directive[Thick],GridLines->{None,{0}},ImageSize->Large]],Column[{"spline","interpolation"},Center]->Show[ListPlot[data,PlotStyle->Directive[Red],PlotMarkers->{Automatic,Large},GridLines->{None,{0}},ImageSize->Large],Plot[sfun[x],{x,0,6},PlotStyle->Directive[Thick],GridLines->{None,{0}},ImageSize->Large]],Column[{"fit","cubic","ploynomial"},Center]->Show[ListPlot[data,PlotStyle->Directive[Red],PlotMarkers->{Automatic,Large},GridLines->{None,{0}},ImageSize->Large],Plot[Evaluate[LinearModelFit[data,{1,x,(x^2),(x^3)},x][x]],{x,0,6},PlotStyle->Directive[Thick],GridLines->{None,{0}},ImageSize->Large]]
}]
]

example[2]:=DynamicModule[{data={{0,0},{1,0},{2,0},{3,0},{4,0},{5,1},{6,10},{7,80},{8,100},{9,150}},hfun,sfun,legend,n},
legend=LineLegend[{ColorData[97][1],ColorData[97][2]},{"Hermite","Spline"}];
Manipulate[
hfun=Interpolation[data,InterpolationOrder->n];
sfun=Interpolation[data,Method->"Spline"];

Show[ListPlot[data,PlotStyle->Directive[Red],PlotMarkers->{Automatic,Medium},GridLines->{None,{0}},ImageSize->Large,PlotRangePadding->Scaled[0.1],PlotRange->All],Plot[{(*ifun[x],*)hfun[x],sfun[x]},{x,0,9},PlotStyle->Directive[Thick],GridLines->{None,{0}},ImageSize->Large,PlotRange->All],
Epilog->{Inset[Show[ListPlot[data,PlotStyle->Directive[Red],PlotMarkers->{Automatic,Medium},GridLines->{None,{0}},ImageSize->200,PlotRange->{{3.5,6.5},5{-1,3}},AspectRatio->1.0],Plot[{(*ifun[x],*)hfun[x],sfun[x]},{x,0,9},PlotStyle->Directive[Thick],GridLines->{None,{0}}]],{0,50},{3,-1}]}
]//Legended[#,legend]&,

{{n,1,"order"},Range[5],ImageSize->Large}
],
TrackedSymbols->{n}
]

example[3]:=Module[{data={{0,0},{1,0},{2,0},{3,0},{4,0},{5,1},{6,10},{7,80},{8,100},{9,150}},hfun,sfun,afun,legend,n=3},
legend=LineLegend[{ColorData[97][1],ColorData[97][2],ColorData[97][3]},{"Hermite","Spline","Akima"}];
hfun=Interpolation[data,InterpolationOrder->n];
sfun=Interpolation[data,Method->"Spline"];
afun=ResourceFunction["AkimaInterpolation"][data];

Show[ListPlot[data,PlotStyle->Directive[Red],PlotMarkers->{Automatic,Medium},GridLines->{None,{0}},ImageSize->Large,PlotRangePadding->Scaled[0.1],PlotRange->All],Plot[{hfun[x],sfun[x],afun[x]},{x,0,9},PlotStyle->Directive[Thick],GridLines->{None,{0}},ImageSize->Large,PlotRange->All],
Epilog->{Inset[Show[ListPlot[data,PlotStyle->Directive[Red],PlotMarkers->{Automatic,Medium},GridLines->{None,{0}},ImageSize->200,PlotRange->{{3.5,6.5},5{-1,3}},AspectRatio->1.0],Plot[{hfun[x],sfun[x],afun[x]},{x,0,9},PlotStyle->Directive[Thick],GridLines->{None,{0}}]],{0,50},{3,-1}]}
]//Legended[#,legend]&
]

theCode[]:=Pane[Style[,FontSize->14],Scrollbars->{False, Automatic},ImageSize->72{10 ,6}]

example[4]:=Join[Table[Grid[data,Frame->{None,None,{{1,i}->True,{1,i+2}->True}}],{i,1,3}],Table[Grid[data,Frame->{None,None,{{i,1}->True,{i+2,1}->True}}],{i,1,3}]]//FlipView