See Fastest data structure for a priority queue - General Usage / Performance - JuliaLang
Interestingly, the break-even for binary search appears close to the MWE’s vector length 20 for xp and yp:
n = 20
271.237 μs (2 allocations: 78.23 KiB)
274.209 μs (2 allocations: 78.23 KiB)
n = 200
837.772 μs (2 allocations: 78.23 KiB)
356.019 μs (2 allocations: 78.23 KiB)
n = 2000
4.551 ms (2 allocations: 78.23 KiB)
523.962 μs (2 allocations: 78.23 KiB)
function weno4_opt(xs, xp, fp)
Ngrid = size(xp)[1]
ε = 1e-6
fs = similar(xs)
left = fp[1]
right = fp[end]
prevβidx = -1
β2 = 0.0
β3 = 0.0
@inbounds for (idx, x) in enumerate(xs)
if x < xp[1]
fs[idx] = left
continue
elseif x > xp[end]
fs[idx] = right
continue
else
i = 2
while @inbounds x > xp[i]
i += 1
end
i -= 1
end
if i == Ngrid
i -= 1
end
if i == 1
xi = xp[i]
xip = xp[i+1]
xipp = xp[i+2]
hi = xip - xi
hip = xipp - xip
yi = fp[i]
yip = fp[i+1]
yipp = fp[i+2]
q3 = yi * ((x - xip) * (x - xipp)) / (hi * (hi + hip))
q3 -= yip * ((x - xi) * (x - xipp)) / (hi * hip)
q3 += yipp * ((x - xi) * (x - xip)) / ((hi + hip) * hip)
fs[idx] = q3
continue
elseif i == Ngrid - 1
xim = xp[i-1]
xi = xp[i]
xip = xp[i+1]
him = xi - xim
hi = xip - xi
yim = fp[i-1]
yi = fp[i]
yip = fp[i+1]
q2 = yim * ((x - xi) * (x - xip)) / (him * (him + hi))
q2 -= yi * ((x - xim) * (x - xip)) / (him * hi)
q2 += yip * ((x - xim) * (x - xi)) / ((him + hi) * hi)
fs[idx] = q2
continue
end
xim = xp[i-1]
xi = xp[i]
xip = xp[i+1]
xipp = xp[i+2]
him = xi - xim
hi = xip - xi
hip = xipp - xip
yim = fp[i-1]
yi = fp[i]
yip = fp[i+1]
yipp = fp[i+2]
q2 = yim * ((x - xi) * (x - xip)) / (him * (him + hi))
q2 -= yi * ((x - xim) * (x - xip)) / (him * hi)
q2 += yip * ((x - xim) * (x - xi)) / ((him + hi) * hi)
q3 = yi * ((x - xip) * (x - xipp)) / (hi * (hi + hip))
q3 -= yip * ((x - xi) * (x - xipp)) / (hi * hip)
q3 += yipp * ((x - xi) * (x - xip)) / ((hi + hip) * hip)
if i != prevβidx
H = him + hi + hip
yyim = - ((2*him + hi)*H + him*(him + hi)) / (him*(him + hi)*H) * yim
yyim += ((him + hi)*H) / (him*hi*(hi + hip)) * yi
yyim -= (him*H) / ((him + hi)*hi*hip) * yip
yyim += (him*(him + hi)) / ((hi + hip)*hip*H) * yipp
yyi = - (hi*(hi + hip)) / (him*(him + hi)*H) * yim
yyi += (hi*(hi + hip) - him*(2*hi + hip)) / (him*hi*(hi + hip)) * yi
yyi += (him*(hi + hip)) / ((him + hi)*hi*hip) * yip
yyi -= (him*hi) / ((hi + hip)*hip*H) * yipp
yyip = (hi*hip) / (him*(him + hi)*H) * yim
yyip -= (hip*(him + hi)) / (him*hi*(hi + hip)) * yi
yyip += ((him + 2*hi)*hip - (him + hi)*hi) / ((him + hi)*hi*hip) * yip
yyip += ((him + hi)*hi) / ((hi + hip)*hip*H) * yipp
yyipp = - ((hi + hip)*hip) / (him*(him + hi)*H) * yim
yyipp += (hip*H) / (him*hi*(hi + hip)) * yi
yyipp -= ((hi + hip) * H) / ((him + hi)*hi*hip) * yip
yyipp += ((2*hip + hi)*H + hip*(hi + hip)) / ((hi + hip)*hip*H) * yipp
β2 = (hi + hip)^2 * (abs(yyip - yyi) / hi - abs(yyi - yyim) / him)^2
β3 = (him + hi)^2 * (abs(yyipp - yyip) / hip - abs(yyip - yyi) / hi)^2
prevβidx = i
end
γ2 = - (x - xipp) / (xipp - xim)
γ3 = (x - xim) / (xipp - xim)
α2 = γ2 / (ε + β2)
α3 = γ3 / (ε + β3)
ω2 = α2 / (α2 + α3)
ω3 = α3 / (α2 + α3)
fs[idx] = ω2*q2 + ω3*q3
end
return fs
end
function weno4_opt_binsearch(xs, xp, fp)
Ngrid = size(xp)[1]
ε = 1e-6
fs = similar(xs)
left = fp[1]
right = fp[end]
prevβidx = -1
β2 = 0.0
β3 = 0.0
@inbounds for (idx, x) in enumerate(xs)
if x < xp[1]
fs[idx] = left
continue
elseif x > xp[end]
fs[idx] = right
continue
else
### linear search
#i = 2
#while @inbounds x > xp[i]
# i += 1
#end
#i -= 1
### binary search
ileft = 2
i = Ngrid - 1
while ileft <= i
mid = (ileft + i) >> 1
if @inbounds x < xp[mid]
i = mid - 1
else
ileft = mid + 1
end
end
end
if i == Ngrid
i -= 1
end
if i == 1
xi = xp[i]
xip = xp[i+1]
xipp = xp[i+2]
hi = xip - xi
hip = xipp - xip
yi = fp[i]
yip = fp[i+1]
yipp = fp[i+2]
q3 = yi * ((x - xip) * (x - xipp)) / (hi * (hi + hip))
q3 -= yip * ((x - xi) * (x - xipp)) / (hi * hip)
q3 += yipp * ((x - xi) * (x - xip)) / ((hi + hip) * hip)
fs[idx] = q3
continue
elseif i == Ngrid - 1
xim = xp[i-1]
xi = xp[i]
xip = xp[i+1]
him = xi - xim
hi = xip - xi
yim = fp[i-1]
yi = fp[i]
yip = fp[i+1]
q2 = yim * ((x - xi) * (x - xip)) / (him * (him + hi))
q2 -= yi * ((x - xim) * (x - xip)) / (him * hi)
q2 += yip * ((x - xim) * (x - xi)) / ((him + hi) * hi)
fs[idx] = q2
continue
end
xim = xp[i-1]
xi = xp[i]
xip = xp[i+1]
xipp = xp[i+2]
him = xi - xim
hi = xip - xi
hip = xipp - xip
yim = fp[i-1]
yi = fp[i]
yip = fp[i+1]
yipp = fp[i+2]
q2 = yim * ((x - xi) * (x - xip)) / (him * (him + hi))
q2 -= yi * ((x - xim) * (x - xip)) / (him * hi)
q2 += yip * ((x - xim) * (x - xi)) / ((him + hi) * hi)
q3 = yi * ((x - xip) * (x - xipp)) / (hi * (hi + hip))
q3 -= yip * ((x - xi) * (x - xipp)) / (hi * hip)
q3 += yipp * ((x - xi) * (x - xip)) / ((hi + hip) * hip)
if i != prevβidx
H = him + hi + hip
yyim = - ((2*him + hi)*H + him*(him + hi)) / (him*(him + hi)*H) * yim
yyim += ((him + hi)*H) / (him*hi*(hi + hip)) * yi
yyim -= (him*H) / ((him + hi)*hi*hip) * yip
yyim += (him*(him + hi)) / ((hi + hip)*hip*H) * yipp
yyi = - (hi*(hi + hip)) / (him*(him + hi)*H) * yim
yyi += (hi*(hi + hip) - him*(2*hi + hip)) / (him*hi*(hi + hip)) * yi
yyi += (him*(hi + hip)) / ((him + hi)*hi*hip) * yip
yyi -= (him*hi) / ((hi + hip)*hip*H) * yipp
yyip = (hi*hip) / (him*(him + hi)*H) * yim
yyip -= (hip*(him + hi)) / (him*hi*(hi + hip)) * yi
yyip += ((him + 2*hi)*hip - (him + hi)*hi) / ((him + hi)*hi*hip) * yip
yyip += ((him + hi)*hi) / ((hi + hip)*hip*H) * yipp
yyipp = - ((hi + hip)*hip) / (him*(him + hi)*H) * yim
yyipp += (hip*H) / (him*hi*(hi + hip)) * yi
yyipp -= ((hi + hip) * H) / ((him + hi)*hi*hip) * yip
yyipp += ((2*hip + hi)*H + hip*(hi + hip)) / ((hi + hip)*hip*H) * yipp
β2 = (hi + hip)^2 * (abs(yyip - yyi) / hi - abs(yyi - yyim) / him)^2
β3 = (him + hi)^2 * (abs(yyipp - yyip) / hip - abs(yyip - yyi) / hi)^2
prevβidx = i
end
γ2 = - (x - xipp) / (xipp - xim)
γ3 = (x - xim) / (xipp - xim)
α2 = γ2 / (ε + β2)
α3 = γ3 / (ε + β3)
ω2 = α2 / (α2 + α3)
ω3 = α3 / (α2 + α3)
fs[idx] = ω2*q2 + ω3*q3
end
return fs
end
using BenchmarkTools
n = 20
xp = sort(randn(n)) # mimicking the OP's data
yp = -sort(-rand(n))
x = collect(LinRange(minimum(xp), maximum(xp), 10001))
@btime weno4_opt(x, xp, yp) evals=1000
@btime weno4_opt_binsearch(x, xp, yp) evals=1000