Skip to content

Commit

Permalink
Merge pull request #66 from gerlero/format
Browse files Browse the repository at this point in the history
Format code with JuliaFormatter
  • Loading branch information
gerlero authored Dec 20, 2023
2 parents 2b8a250 + 5846dad commit 7cb27a1
Show file tree
Hide file tree
Showing 6 changed files with 121 additions and 70 deletions.
1 change: 1 addition & 0 deletions .JuliaFormatter
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
style = "sciml"
8 changes: 8 additions & 0 deletions .github/workflows/Format.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
name: Format suggestions
on:
pull_request:
jobs:
code-style:
runs-on: ubuntu-latest
steps:
- uses: julia-actions/julia-format@v2
1 change: 1 addition & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
[![Build Status](https://github.com/gerlero/PCHIPInterpolation.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/gerlero/PCHIPInterpolation.jl/actions/workflows/CI.yml?query=branch%3Amain)
[![Coverage](https://codecov.io/gh/gerlero/PCHIPInterpolation.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/gerlero/PCHIPInterpolation.jl)
[![PkgEval](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/P/PCHIPInterpolation.svg)](https://JuliaCI.github.io/NanosoldierReports/pkgeval_badges/P/PCHIPInterpolation.html)
[![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)

PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) spline interpolation of arbitrarily spaced one-dimensional data in Julia. This package is a fork of [SimplePCHIP](https://github.com/slabanja/SimplePCHIP) with some extra features.

Expand Down
6 changes: 3 additions & 3 deletions ext/PCHIPInterpolationRecipesBaseExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,11 @@ module PCHIPInterpolationRecipesBaseExt
using PCHIPInterpolation
using RecipesBase

@recipe function _(itp::Interpolator; markershape=:none) \
@recipe function _(itp::Interpolator; markershape = :none)
@series begin
markershape := :none
plotdensity = clamp(10*length(itp.xs), 1000, 100000)
x = range(first(itp.xs), last(itp.xs), length=plotdensity)
plotdensity = clamp(10 * length(itp.xs), 1000, 100000)
x = range(first(itp.xs), last(itp.xs), length = plotdensity)
return x, itp.(x)
end
if markershape !== :none
Expand Down
55 changes: 36 additions & 19 deletions src/PCHIPInterpolation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,12 @@ module PCHIPInterpolation
export Interpolator, integrate


_is_strictly_increasing(xs::AbstractVector) = all(a < b for (a,b) in zip(xs[begin:end-1], xs[begin+1:end]))
_is_strictly_increasing(xs::AbstractVector) =
all(a < b for (a, b) in zip(xs[begin:end-1], xs[begin+1:end]))
_is_strictly_increasing(xs::AbstractRange) = step(xs) > 0

function _pchip_edge_derivative(h1, h2, Δ1, Δ2)
d = ((2h1 + h2)*Δ1 - h2*Δ2)/(h1 + h2)
d = ((2h1 + h2) * Δ1 - h2 * Δ2) / (h1 + h2)
if sign(d) != sign(Δ1)
d = zero(d)
elseif sign(Δ1) != sign(Δ2) && abs(d) > abs(3Δ1)
Expand All @@ -18,9 +19,9 @@ end

function _pchip_ds_scipy(xs::AbstractVector, ys::AbstractVector)
h(i) = xs[i+1] - xs[i]
Δ(i) = (ys[i+1] - ys[i])/h(i)
Δ(i) = (ys[i+1] - ys[i]) / h(i)

ds = similar(ys./xs)
ds = similar(ys ./ xs)

is = eachindex(xs, ys, ds)

Expand All @@ -37,13 +38,23 @@ function _pchip_ds_scipy(xs::AbstractVector, ys::AbstractVector)
else
wl = 2hr + hl
wr = hr + 2hl
ds[i] = (wl + wr)/(wl/Δl + wr/Δr)
ds[i] = (wl + wr) / (wl / Δl + wr / Δr)
end
Δl = Δr
hl = hr
end
ds[begin] = _pchip_edge_derivative(h(first(is)), h(first(is)+1), Δ(first(is)), Δ(first(is)+1))
ds[end] = _pchip_edge_derivative(h(last(is)-1), h(last(is)-2), Δ(last(is)-1), Δ(last(is)-2))
ds[begin] = _pchip_edge_derivative(
h(first(is)),
h(first(is) + 1),
Δ(first(is)),
Δ(first(is) + 1),
)
ds[end] = _pchip_edge_derivative(
h(last(is) - 1),
h(last(is) - 2),
Δ(last(is) - 1),
Δ(last(is) - 2),
)
end
return ds
end
Expand All @@ -54,17 +65,21 @@ struct Interpolator{Xs,Ys,Ds}
ds::Ds

function Interpolator(xs::AbstractVector, ys::AbstractVector)
length(eachindex(xs, ys)) 2 || throw(ArgumentError("inputs must have at least 2 elements"))
_is_strictly_increasing(xs) || throw(ArgumentError("xs must be strictly increasing"))
length(eachindex(xs, ys)) 2 ||
throw(ArgumentError("inputs must have at least 2 elements"))
_is_strictly_increasing(xs) ||
throw(ArgumentError("xs must be strictly increasing"))

ds = _pchip_ds_scipy(xs, ys)

new{typeof(xs),typeof(ys),typeof(ds)}(xs, ys, ds)
end

function Interpolator(xs::AbstractVector, ys::AbstractVector, ds::AbstractVector)
length(eachindex(xs, ys, ds)) 2 || throw(ArgumentError("inputs must have at least 2 elements"))
_is_strictly_increasing(xs) || throw(ArgumentError("xs must be strictly increasing"))
length(eachindex(xs, ys, ds)) 2 ||
throw(ArgumentError("inputs must have at least 2 elements"))
_is_strictly_increasing(xs) ||
throw(ArgumentError("xs must be strictly increasing"))

new{typeof(xs),typeof(ys),typeof(ds)}(xs, ys, ds)
end
Expand Down Expand Up @@ -103,7 +118,7 @@ end
return imax - 1 # Treat right endpoint as part of rightmost interval
end

i = imin + (imax - imin + 1)÷2
i = imin + (imax - imin + 1) ÷ 2
while imin < imax
if x < @inbounds xs[i]
imax = i - 1
Expand All @@ -112,7 +127,7 @@ end
else
break
end
i = imin + (imax - imin + 1)÷2
i = imin + (imax - imin + 1) ÷ 2
end

return i
Expand Down Expand Up @@ -177,10 +192,10 @@ function _evaluate(itp::Interpolator, x, i)
d1 = _derivative(itp, Val(:begin), i)
d2 = _derivative(itp, Val(:end), i)

return (y1 * ((x2-x)/h)
+ y2 * ((x-x1)/h)
- d1*h * ((x2-x)/h)
+ d2*h * ((x-x1)/h))
return (
y1 * ((x2 - x) / h) + y2 * ((x - x1) / h) - d1 * h * ((x2 - x) / h) +
d2 * h * ((x - x1) / h)
)
end

@inline _evaluate(itp::Interpolator, x) = _evaluate(itp, x, _findinterval(itp, x))
Expand All @@ -191,7 +206,9 @@ end
@inline function _integrate(itp::Interpolator, a, b, i)
a_ = _x(itp, a, i)
b_ = _x(itp, b, i)
return (b_ - a_)/6*(_evaluate(itp, a, i) + 4*_evaluate(itp, (a_ + b_)/2, i) + _evaluate(itp, b, i)) # Simpson's rule
return (b_ - a_) / 6 * (
_evaluate(itp, a, i) + 4 * _evaluate(itp, (a_ + b_) / 2, i) + _evaluate(itp, b, i)
) # Simpson's rule
end

@inline function _integrate(itp::Interpolator, a, b, i, j)
Expand All @@ -200,7 +217,7 @@ end
end

integral = _integrate(itp, a, Val(:end), i)
for k in i+1:j-1
for k = i+1:j-1
integral += _integrate(itp, Val(:begin), Val(:end), k)
end
integral += _integrate(itp, Val(:begin), b, j)
Expand Down
Loading

0 comments on commit 7cb27a1

Please sign in to comment.