Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
######################################################################
using PyPlot
using Grid
using Dierckx
using NumericalMath
using splines
#using ApproXD
meshgrid(x::Vector, y::Vector) = (repmat(x, 1, length(y))',
repmat(y, 1, length(x)))
######################################################################
# Function to interpolate
f(x::Float64) = -(x^(-1.0))
# Grids for interpolation (regular/irregular)
zmin = -0.1
zmax = 1000.0
ymin = 0.2
ymax = 30.0
xmin = ymin + zmin
xmax = ymax + zmax
xpoints = 100
ypoints = 20
zpoints = 15
xgrid = linspace(xmin, xmax, xpoints)
ygrid = linspace(ymin, ymax, ypoints)
zgrid = linspace(zmin, zmax, zpoints)
power = 3 # Curvature of the irregular grids
x_irrgrid = zeros(xpoints)
x_irrgridexp = zeros(xpoints)
xdistexp = (xmax - xmin)^(1/power)
xinc = xdistexp/(xpoints-1)
for i = 1: xpoints
x_irrgridexp[i] = (i-1)*xinc
end
x_irrgrid[:] = x_irrgridexp.^power + xmin
y_irrgrid = zeros(ypoints)
y_irrgridexp = zeros(ypoints)
ydisteyp = (ymax - ymin)^(1/power)
yinc = ydisteyp/(ypoints-1)
for i = 1: ypoints
y_irrgridexp[i] = (i-1)*yinc
end
y_irrgrid[:] = y_irrgridexp.^power + ymin
z_irrgrid = zeros(zpoints)
z_irrgridexp = zeros(zpoints)
zdistezp = (zmax - zmin)^(1/power)
zinc = zdistezp/(zpoints-1)
for i = 1: zpoints
z_irrgridexp[i] = (i-1)*zinc
end
z_irrgrid[:] = z_irrgridexp.^power + zmin
######################################################################
# Evaluate functions on grid points
V_1D = zeros(xpoints)
V_1D_irr = zeros(xpoints)
V_2D = zeros(ypoints, zpoints)
V_2D_irr = zeros(ypoints, zpoints)
# Values on grid, 1 dimension
for x = 1:xpoints
xt = xgrid[x]
xt2 = x_irrgrid[x]
V_1D[x] = f(xt)
V_1D_irr[x] = f(xt2)
end
# Values on grid, 2 dimensions
for z = 1:zpoints
for y = 1:ypoints
zt = zgrid[z]
yt = ygrid[y]
zt2 = z_irrgrid[z]
yt2 = y_irrgrid[y]
V_2D[y, z, 1] = f(yt+zt)
V_2D_irr[y, z, 1] = f(yt2 + zt2)
end
end
######################################################################
# Interpolate in 1 dimension
xrange = range(xgrid[1], xgrid[2]-xgrid[1], xpoints)
# Grid
V_1D_int_Lin = CoordInterpGrid(xrange, V_1D, BCnearest, InterpLinear)
V_1D_int_Quad = CoordInterpGrid(xrange, V_1D, BCnearest, InterpQuadratic)
#V_1D_int_Cub = CoordInterpGrid(xrange, V_1D[:, 1], BCnil, InterpCubic)
# ApproXD
# V_1D_int_ApproXD =
# Dierckx (regular/irregular grid)
V_1D_int_Dierckx = Spline1D(xgrid, V_1D)
V_1D_int_irr_Dierckx = Spline1D(x_irrgrid, V_1D_irr)
######################################################################
# Interpolate in 2 dimensions
yrange = range(ygrid[1], ygrid[2]-ygrid[1], ypoints)
zrange = range(zgrid[1], zgrid[2]-zgrid[1], zpoints)
# Grid
V_2D_int_Lin = CoordInterpGrid((yrange, zrange), V_2D, BCnearest, InterpLinear)
V_2D_int_Quad = CoordInterpGrid((yrange, zrange), V_2D, BCnearest, InterpQuadratic)
#V_2D_int_Cub = CoordInterpGrid((yrange, zrange), V_2D, BCnan, InterpCubic)
# ApproXD
# V_2D_int_ApproXD =
# Dierckx (regular/irregular grid)
V_2D_int_Dierckx = Spline2D(ygrid, zgrid, V_2D)
V_2D_int_irr_Dierckx = Spline2D(y_irrgrid, z_irrgrid, V_2D_irr)
# Econforge/splines
CubicSplines = interpolant_cspline([ymin,zmin], [ymax,zmax], [ypoints,zpoints], V_2D )
# NumericalMath - Piecewise cubic hermitean interpolating Polynomial
# As in Matlab, NumericalMath doesn't separate interpolation and
# evaluation, there is only one call to pchip(xgrid, V_1D, x)
######################################################################
# Evaluate Regular and Irregular Grid Interpolants Off-Grid
xop = 5*xpoints # Number of off-grid points
yop = 5*ypoints
zop = 5*zpoints
xoffgrid = linspace(xgrid[1], xgrid[end], xop)
yoffgrid = linspace(ygrid[1], ygrid[end], yop)
zoffgrid = linspace(zgrid[1], zgrid[end], zop)
V_1D_offgrid_Lin = zeros(xop)
V_1D_offgrid_Quad = zeros(xop)
V_1D_offgrid_Dierckx = zeros(xop)
V_1D_offgrid_pchip = zeros(xop)
V_1D_irr_offgrid_Dierckx = zeros(xop)
V_1D_actual = zeros(xop)
V_2D_offgrid_Lin = zeros(yop, zop)
V_2D_offgrid_Quad = zeros(yop, zop)
V_2D_offgrid_Dierckx = zeros(yop, zop)
V_2D_offgrid_CSplines = zeros(yop, zop)
V_2D_irr_offgrid_Dierckx = zeros(yop, zop)
V_2D_actual = zeros(yop, zop)
for x = 1:xop
xt = xoffgrid[x]
# Grid
V_1D_offgrid_Lin[x] = V_1D_int_Lin[xt]
V_1D_offgrid_Quad[x] = V_1D_int_Quad[xt]
# Dierckx (regular/irregular grid)
V_1D_offgrid_Dierckx[x] = evaluate(V_1D_int_Dierckx, xt)
V_1D_irr_offgrid_Dierckx[x] = evaluate(V_1D_int_irr_Dierckx, xt)
# NumericalMath
V_1D_offgrid_pchip[x] = pchip(xgrid, V_1D, xt)
# Actual Function value
V_1D_actual[x] = f(xt)
end
for y = 1:yop
for z = 1:zop
yt = yoffgrid[y]
zt = zoffgrid[z]
# Grid
V_2D_offgrid_Lin[y, z] = V_2D_int_Lin[yt, zt]
V_2D_offgrid_Quad[y, z] = V_2D_int_Quad[yt, zt]
V_2D_offgrid_CSplines[y, z] = CubicSplines(yt, zt)[1]
# Dierckx (regular/irregular grid)
V_2D_offgrid_Dierckx[y, z] = evaluate(V_2D_int_Dierckx, yt, zt)
V_2D_irr_offgrid_Dierckx[y, z] = evaluate(V_2D_int_irr_Dierckx, yt, zt)
# Actual Function value
V_2D_actual[y, z] = f(yt + zt)
end
end
######################################################################
# Calculate Errors
Err_1D_Lin = (V_1D_actual - V_1D_offgrid_Lin)./V_1D_actual
Err_1D_Quad = (V_1D_actual - V_1D_offgrid_Quad)./V_1D_actual
Err_1D_Dierckx = (V_1D_actual - V_1D_offgrid_Dierckx)./V_1D_actual
Err_1D_pchip = (V_1D_actual - V_1D_offgrid_pchip)./V_1D_actual
Err_1D_irr_Dierckx = (V_1D_actual - V_1D_irr_offgrid_Dierckx)./V_1D_actual
Err_2D_Lin = (V_2D_actual - V_2D_offgrid_Lin)./V_2D_actual
Err_2D_Quad = (V_2D_actual - V_2D_offgrid_Quad)./V_2D_actual
Err_2D_CSplines = (V_2D_actual - V_2D_offgrid_CSplines)./V_2D_actual
Err_2D_Dierckx = (V_2D_actual - V_2D_offgrid_Dierckx)./V_2D_actual
Err_2D_irr_Dierckx = (V_2D_actual - V_2D_irr_offgrid_Dierckx)./V_2D_actual
@printf "∑(abs(ɛ)) for 1D linear interpolation with Grid is %.3f\n" sum(abs(Err_1D_Lin))/xop
@printf "∑(abs(ɛ)) for 1D quadratic interpolation with Grid is %.3f\n" sum(abs(Err_1D_Quad))/xop
@printf "∑(abs(ɛ)) for 1D spline interpolation on a regular grid with Dierckx is %.3f\n" sum(abs(Err_1D_Dierckx))/xop
@printf "∑(abs(ɛ)) for 1D spline interpolation on an irregular grid with Dierckx is %.3f\n" sum(abs(Err_1D_irr_Dierckx))/xop
@printf "∑(abs(ɛ)) for 1D pchip interpolation with NumericalMath is %.2f\n\n" sum(abs(Err_1D_pchip))/xop
@printf "∑(abs(ɛ)) for 2D linear interpolation with Grid is %.3f\n" sum(abs(Err_2D_Lin))/(yop*zop)
@printf "∑(abs(ɛ)) for 2D quadratic interpolation with Grid is %.3f\n" sum(abs(Err_2D_Quad))/(yop*zop)
@printf "∑(abs(ɛ)) for 2D cubic interpolation with Cubic Splines is %.3f\n" sum(abs(Err_2D_CSplines))/(yop*zop)
@printf "∑(abs(ɛ)) for 2D spline interpolation on a regular grid with Dierckx is %.3f\n" sum(abs(Err_2D_Dierckx))/(yop*zop)
@printf "∑(abs(ɛ)) for 2D spline interpolation on an irregular grid with Dierckx is %.3f\n" sum(abs(Err_2D_irr_Dierckx))/(yop*zop)
######################################################################
# Plots
fig, ax = PyPlot.subplots(2,3)
ax[1,1][:plot](V_1D_actual[1:20], label = "True Value")
ax[1,1][:plot](V_1D_offgrid_Lin[1:20], label = "Grid Linearly Interpolated Value")
ax[1,1][:plot](Err_1D_Lin[1:20], label = "Interpolation Error")
ax[1,1][:legend](loc = 0)
ax[2,1][:plot](V_1D_actual[1:20], label = "True Value")
ax[2,1][:plot](V_1D_offgrid_Quad[1:20], label = "Grid Quadratically Interpolated Value")
ax[2,1][:plot](Err_1D_Quad[1:20], label = "Interpolation Error")
ax[2,1][:legend](loc = 0)
ax[1,2][:plot](V_1D_actual[1:20], label = "True Value")
ax[1,2][:plot](V_1D_offgrid_Dierckx[1:20], label = "Dierckx Interpolated Value (regular grid)")
ax[1,2][:plot](Err_1D_Dierckx[1:20], label = "Interpolation Error")
ax[1,2][:legend](loc = 0)
ax[2,2][:plot](V_1D_actual[1:20], label = "True Value")
ax[2,2][:plot](V_1D_offgrid_pchip[1:20], label = "PCHIP Interpolated Value")
ax[2,2][:plot](Err_1D_pchip[1:20], label = "Interpolation Error")
ax[2,2][:legend](loc = 0)
ax[2,3][:plot](V_1D_actual[1:20], label = "True Value")
ax[2,3][:plot](V_1D_irr_offgrid_Dierckx[1:20], label = "Dierckx Interpolated Value (irregular grid)")
ax[2,3][:plot](Err_1D_irr_Dierckx[1:20], label = "Interpolation Error")
ax[2,3][:legend](loc = 0)
fig[:suptitle]("Interpolation in 1D")
plt.show()
fig2 = figure()
ax = fig2[:add_subplot](231, projection="3d")
yg, zg = meshgrid(yoffgrid[1:ypoints], zoffgrid[1:zpoints])
vg = V_2D_offgrid_Lin[1:ypoints, 1:zpoints]'
ax[:plot_surface](yg, zg, vg, rstride = 1, cstride = 1,
cmap=ColorMap("jet"), alpha=0.5, linewidth=0.25)
ax[:set_title]("Grid Linear Interpolation")
ax = fig2[:add_subplot](232, projection="3d")
yg, zg = meshgrid(yoffgrid[1:ypoints], zoffgrid[1:zpoints])
vg = V_2D_offgrid_Quad[1:ypoints, 1:zpoints]'
ax[:plot_surface](yg, zg, vg, rstride = 1, cstride = 1,
cmap=ColorMap("jet"), alpha=0.5, linewidth=0.25)
ax[:set_title]("Grid Quadratic Interpolation")
ax = fig2[:add_subplot](233, projection="3d")
yg, zg = meshgrid(yoffgrid[1:ypoints], zoffgrid[1:zpoints])
vg = V_2D_offgrid_CSplines[1:ypoints, 1:zpoints]'
ax[:plot_surface](yg, zg, vg, rstride = 1, cstride = 1,
cmap=ColorMap("jet"), alpha=0.5, linewidth=0.25)
ax[:set_title]("Grid Cubic Interpolation")
ax = fig2[:add_subplot](234, projection="3d")
yg, zg = meshgrid(yoffgrid[1:ypoints], zoffgrid[1:zpoints])
vg = V_2D_offgrid_Dierckx[1:20, 1:15]'
ax[:plot_surface](yg, zg, vg, rstride = 1, cstride = 1,
cmap=ColorMap("jet"), alpha=0.5, linewidth=0.25)
ax[:set_title]("Dierckx Spline Interpolation (regular grid")
ax = fig2[:add_subplot](235, projection="3d")
yg, zg = meshgrid(yoffgrid[1:ypoints], zoffgrid[1:zpoints])
vg = V_2D_irr_offgrid_Dierckx[1:ypoints, 1:zpoints]'
ax[:plot_surface](yg, zg, vg, rstride = 1, cstride = 1,
cmap=ColorMap("jet"), alpha=0.5, linewidth=0.25)
ax[:set_title]("Dierckx Spline Interpolation (irregular grid)")
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment