Interpolations with Unitful Vectors

General Usage
,
1

Often I need to perform interpolations of unit-aware vectors and I would prefer to have the flexibility to use the interpolation functions offered by Dierckx, Interpolations etc.

I haven’t found this in docs or discourse so I wanted to provide it here in case it’s useful to anybody else and see if maybe there’s already a better way that I overlooked.

using Unitful, Interpolations

# Creating an x and y vector with units
Xs = (0:0.1:5) .* u"m"
Ys = Xs.^2

# Defining a function which will return an interpolator which strips units according to input xs and applies units to output interpolated values based on the units of ys
function UnitfulLinterp(xs, ys; interpolator_method=LinearInterpolation)
	x_unit = unit(xs[1])
	y_unit = unit(ys[1])

	xs = ustrip.(x_unit, xs)
	ys = ustrip.(y_unit, ys)

	itp = interpolator_method(xs, ys)

	interpolator(x) = itp(ustrip(x_unit, x)) * y_unit
	
	return interpolator
end

# creating my unit-aware interpolator, itp
itp = UnitfulLinterp(Xs, Ys)

# passing a single value to be interpolated
x_to_eval = 5u"mm"
interpolated_y = itp(x_to_eval)
> 0.0005008241058144431 m^2

# passing a vector of values and broadcasting the interpolator
xs_to_eval = [5.0u"mm", 12.0u"inch", 15u"ft"]
interpolated_ys = itp.(xs_to_eval)
> 0.000500824 m^2
> 0.0933604 m^2
> 20.9052 m^2
2 Likes
2

I think this is an oversight in Interpolations. Unitful interpolation works fine for Gridded interpolation as shown in the tests, the trick there is that both positions and values are arrays, not ranges. For your example here:

Xs = collect((0:0.1:5)u"m")
Ys = Xs .^ 2
itp = LinearInterpolation(Xs, Ys)

itp(5.0u"inch")

The problem is that currently scaled interpolation does not work, lets see if this PR to fix it is accepted, then this will work:

Xs = (0:0.1:5)u"m"
Ys = Xs .^ 2
sitp = scale(interpolate(Ys, BSpline(Linear())), Xs)

sitp(5.0u"inch")
2 Likes