From d6346c6925fb2dbedc32dc63c6a4b22cae864b06 Mon Sep 17 00:00:00 2001 From: Alexander Morley Date: Tue, 13 Feb 2018 12:12:54 +0000 Subject: [PATCH 1/3] update --- src/gridded/cubic.jl | 135 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 135 insertions(+) create mode 100644 src/gridded/cubic.jl diff --git a/src/gridded/cubic.jl b/src/gridded/cubic.jl new file mode 100644 index 00000000..321cb9db --- /dev/null +++ b/src/gridded/cubic.jl @@ -0,0 +1,135 @@ +function getindex(itp::GriddedInterpolation{T,N,TCoefs,Interpolations.Gridded{Interpolations.Cubic{Interpolations.Reflect}},K,P}, x::Number) where{T,N,TCoefs,K,P} + a,b,c = coefficients(itp.knots[1], itp.coefs) + interpolate(x, a, b, c, itp.knots[1], itp.coefs) +end + +function getindex(itp::GriddedInterpolation{T,N,TCoefs,Interpolations.Gridded{Interpolations.Cubic{Interpolations.Reflect}},K,P}, x::AbstractVector) where{T,N,TCoefs,K,P} + a,b,c = coefficients(itp.knots[1], itp.coefs) + interpolate.(x, [a], [b], [c], [itp.knots[1]], [itp.coefs]) +end + + +""" + interpolate(x, a, b, c, X, Y, v=false) +Interpoalte at location x using coefficients a,b,c fitted to X & Y. +If i in in the range of x: + Sᵢ(x) = a(x - xᵢ)³ + b(x - xᵢ)² + c(x - xᵢ) + dᵢ +Otherwise extrapolate with a = 0 (i.e. a 1st or 2nd order spline) +""" +function interpolate(x, a, b, c, X, Y, v=false) + idx = max(searchsortedfirst(X,x)-1,1) + h = x - X[idx] # (x - xᵢ) + n = length(X) + + if x < X[1] # extrapolation to the left + return (b[1]*h + c[1])*h + Y[1]; + elseif x > X[n] # extrapolation to the right + return (b[n]*h + c[n])*h + Y[n]; + else # interpolation + return ((a[idx]*h + b[idx])*h + c[idx])*h + Y[idx]; + end + interpol +end + + +function coefficients(x, y; + x1 = 0., + xend = 0., + m_force_linear_extrapolation = true, + boundaries = (:first_deriv, :first_deriv) + ) + n = length(x) + A = zeros(n,n) + rx = zeros(n) + for i in 2:(n-1) + A[i-1,i] = 1.0/3.0*(x[i]-x[i-1]); + A[i,i] = 2.0/3.0*(x[i+1]-x[i-1]); + A[i+1,i] = 1.0/3.0*(x[i+1]-x[i]); + rx[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]) + end + + boundary_start,boundary_end = boundaries + if (boundary_start == :second_deriv) + A[1,1] = 2.0; + A[2,1] = 0.0; + rx[0] = x1; + elseif (boundary_start == :first_deriv) + A[1,1] = 2.0*(x[2]-x[1]); + A[2,1] = 1.0*(x[2]-x[1]); + rx[1] = 3.0*((y[2]-y[1])/(x[2]-x[1])-x1); + else + assert(false); + end + + if (boundary_end == :second_deriv) + A[n,n] = 2.0; + A[n-1,n] = 0.0; + rx[n] = xend; + elseif (boundary_end == :first_deriv) + A[n,n] = 2.0*(x[n]-x[n-1]); + A[n-1,n] =1.0*(x[n]-x[n-1]); + rx[n]=3.0*(xend-(y[n]-y[n-1])/(x[n]-x[n-1])); + end + + b = A \ rx; + + a = zeros(n) + c = zeros(n) + for i in 1:(n-1) + a[i]=1.0/3.0*(b[i+1]-b[i])/(x[i+1]-x[i]) + c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])- 1.0/3.0*(2.0*b[i]+b[i+1])*(x[i+1]-x[i]) + end + + b[1] = (m_force_linear_extrapolation==false) ? b[1] : 0.0; + h = x[n]-x[n-1]; + a[n] = 0.0 + c[n] = 3a[n-1] * (h * h) + 2b[n-1] * h + c[n-1] + m_force_linear_extrapolation && (b[n]=0.0) + return a, b, c +end + +##### + +function define_indices_d(::Type{Gridded{Cubic{Reflect}}}, d, pad) + symix, symixp, symx = Symbol("ix_",d), Symbol("ixp_",d), Symbol("x_",d) + quote + $symix = clamp($symix, 1, size(itp, $d)-1) + $symixp = $symix + 1 + end +end + +function coefficients(::Type{Gridded{Cubic{Reflect}}}, N, d) + symix, symixp, symx = Symbol("ix_",d), Symbol("ixp_",d), Symbol("x_",d) + sym, symp, symfx = Symbol("c_",d), Symbol("cp_",d), Symbol("fx_",d) + symk, symkix = Symbol("k_",d), Symbol("kix_",d) + quote + $symkix = $symk[$symix] + $symfx = ($symx - $symkix)/($symk[$symixp] - $symkix) + $sym = 1 - $symfx + $symp = $symfx + end +end + +function gradient_coefficients(::Type{Gridded{Cubic{Reflect}}}, d) + sym, symp = Symbol("c_",d), Symbol("cp_",d) + symk, symix = Symbol("k_",d), Symbol("ix_",d) + symixp = Symbol("ixp_",d) + quote + $symp = 1/($symk[$symixp] - $symk[$symix]) + $sym = - $symp + end +end + +# This assumes fractional values 0 <= fx_d <= 1, integral values ix_d and ixp_d (typically ixp_d = ix_d+1, +#except at boundaries), and an array itp.coefs +function index_gen(::Type{Gridded{Cubic{Reflect}}}, ::Type{IT}, N::Integer, offsets...) where IT<:DimSpec{Gridded} + if length(offsets) < N + d = length(offsets)+1 + sym = Symbol("c_", d) + symp = Symbol("cp_", d) + return :($sym * $(index_gen(IT, N, offsets..., 0)) + $symp * $(index_gen(IT, N, offsets..., 1))) + else + indices = [offsetsym(offsets[d], d) for d = 1:N] + return :(itp.coefs[$(indices...)]) + end +end From e83e657cfc079c8c64fcc54f85ab87da2320b1cc Mon Sep 17 00:00:00 2001 From: Alexander Morley Date: Mon, 19 Feb 2018 13:32:33 +0000 Subject: [PATCH 2/3] save --- src/gridded/cubic.jl | 82 ++++++++++++++++++++---------------------- src/gridded/gridded.jl | 1 + 2 files changed, 40 insertions(+), 43 deletions(-) diff --git a/src/gridded/cubic.jl b/src/gridded/cubic.jl index 321cb9db..14a98530 100644 --- a/src/gridded/cubic.jl +++ b/src/gridded/cubic.jl @@ -1,43 +1,39 @@ function getindex(itp::GriddedInterpolation{T,N,TCoefs,Interpolations.Gridded{Interpolations.Cubic{Interpolations.Reflect}},K,P}, x::Number) where{T,N,TCoefs,K,P} a,b,c = coefficients(itp.knots[1], itp.coefs) - interpolate(x, a, b, c, itp.knots[1], itp.coefs) + interpolate(x, a, b, c, d, itp.knots[1]) end function getindex(itp::GriddedInterpolation{T,N,TCoefs,Interpolations.Gridded{Interpolations.Cubic{Interpolations.Reflect}},K,P}, x::AbstractVector) where{T,N,TCoefs,K,P} a,b,c = coefficients(itp.knots[1], itp.coefs) - interpolate.(x, [a], [b], [c], [itp.knots[1]], [itp.coefs]) + interpolate.(x, [a], [b], [c], [d], [itp.knots[1]]) end """ - interpolate(x, a, b, c, X, Y, v=false) + interpolate(x, a, b, c, d, X, v=false) Interpoalte at location x using coefficients a,b,c fitted to X & Y. If i in in the range of x: Sᵢ(x) = a(x - xᵢ)³ + b(x - xᵢ)² + c(x - xᵢ) + dᵢ Otherwise extrapolate with a = 0 (i.e. a 1st or 2nd order spline) """ -function interpolate(x, a, b, c, X, Y, v=false) +function interpolate(x, a, b, c, d, X, v=false) idx = max(searchsortedfirst(X,x)-1,1) - h = x - X[idx] # (x - xᵢ) n = length(X) - - if x < X[1] # extrapolation to the left - return (b[1]*h + c[1])*h + Y[1]; - elseif x > X[n] # extrapolation to the right - return (b[n]*h + c[n])*h + Y[n]; - else # interpolation - return ((a[idx]*h + b[idx])*h + c[idx])*h + Y[idx]; + idx = idx >= n ? idx - 1 : idx + if idx < n # interpolation + return a[idx]*(x^3) + b[idx]*(x^2) + c[idx]x + d[idx] end - interpol + error() end function coefficients(x, y; - x1 = 0., - xend = 0., - m_force_linear_extrapolation = true, - boundaries = (:first_deriv, :first_deriv) + force_linear_extrapolation = true, + boundary_condition = :natural ) + const x1 = 0. + const xend = 0. + n = length(x) A = zeros(n,n) rx = zeros(n) @@ -47,28 +43,25 @@ function coefficients(x, y; A[i+1,i] = 1.0/3.0*(x[i+1]-x[i]); rx[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]) end - - boundary_start,boundary_end = boundaries - if (boundary_start == :second_deriv) - A[1,1] = 2.0; - A[2,1] = 0.0; - rx[0] = x1; - elseif (boundary_start == :first_deriv) - A[1,1] = 2.0*(x[2]-x[1]); - A[2,1] = 1.0*(x[2]-x[1]); - rx[1] = 3.0*((y[2]-y[1])/(x[2]-x[1])-x1); + + if boundary_condition == :natural + A[1,1] = 2.0 + A[2,1] = 0.0 + rx[1] = x1 + + A[n,n] = 2.0 + A[n-1,n] = 0.0 + rx[n] = xend + elseif boundary_condition == :periodic + A[1,1] = 2.0*(x[2]-x[1]) + A[2,1] = 1.0*(x[2]-x[1]) + rx[1] = 3.0*((y[2]-y[1])/(x[2]-x[1])-x1) + + A[n,n] = 2.0*(x[n]-x[n-1]) + A[n-1,n] =1.0*(x[n]-x[n-1]) + rx[n]=3.0*(xend-(y[n]-y[n-1])/(x[n]-x[n-1])) else - assert(false); - end - - if (boundary_end == :second_deriv) - A[n,n] = 2.0; - A[n-1,n] = 0.0; - rx[n] = xend; - elseif (boundary_end == :first_deriv) - A[n,n] = 2.0*(x[n]-x[n-1]); - A[n-1,n] =1.0*(x[n]-x[n-1]); - rx[n]=3.0*(xend-(y[n]-y[n-1])/(x[n]-x[n-1])); + error("Boundary condition $boundary_condition not recognised.") end b = A \ rx; @@ -80,11 +73,14 @@ function coefficients(x, y; c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])- 1.0/3.0*(2.0*b[i]+b[i+1])*(x[i+1]-x[i]) end - b[1] = (m_force_linear_extrapolation==false) ? b[1] : 0.0; - h = x[n]-x[n-1]; - a[n] = 0.0 - c[n] = 3a[n-1] * (h * h) + 2b[n-1] * h + c[n-1] - m_force_linear_extrapolation && (b[n]=0.0) + if force_linear_extrapolation + b[1] = 0.0 + b[n] = 0.0 + a[n] = 0.0 + @show a[1] + h = x[n]-x[n-1]; + c[n] = 3a[n-1] * (h * h) + 2b[n-1] * h + c[n-1] + end return a, b, c end diff --git a/src/gridded/gridded.jl b/src/gridded/gridded.jl index 74daab41..9a6bda31 100644 --- a/src/gridded/gridded.jl +++ b/src/gridded/gridded.jl @@ -69,4 +69,5 @@ ubound(itp::GriddedInterpolation, d, inds) = itp.knots[d][end] include("constant.jl") include("linear.jl") +include("cubic.jl") include("indexing.jl") From 6b638ac9ad03ca6a1f7231b362969be843395d88 Mon Sep 17 00:00:00 2001 From: Alexander Morley Date: Mon, 19 Feb 2018 13:37:43 +0000 Subject: [PATCH 3/3] first pass --- src/gridded/cubic.jl | 134 +++++++++++++++++++++++++++---------------- 1 file changed, 85 insertions(+), 49 deletions(-) diff --git a/src/gridded/cubic.jl b/src/gridded/cubic.jl index 14a98530..097ebdc9 100644 --- a/src/gridded/cubic.jl +++ b/src/gridded/cubic.jl @@ -1,10 +1,10 @@ -function getindex(itp::GriddedInterpolation{T,N,TCoefs,Interpolations.Gridded{Interpolations.Cubic{Interpolations.Reflect}},K,P}, x::Number) where{T,N,TCoefs,K,P} - a,b,c = coefficients(itp.knots[1], itp.coefs) +function getindex(itp::GriddedInterpolation{T,N,TCoefs,Interpolations.Gridded{Interpolations.Cubic{Interpolations.Natural}},K,P}, x::Number) where{T,N,TCoefs,K,P} + a,b,c,d = coefficients(itp.knots[1], itp.coefs) interpolate(x, a, b, c, d, itp.knots[1]) end -function getindex(itp::GriddedInterpolation{T,N,TCoefs,Interpolations.Gridded{Interpolations.Cubic{Interpolations.Reflect}},K,P}, x::AbstractVector) where{T,N,TCoefs,K,P} - a,b,c = coefficients(itp.knots[1], itp.coefs) +function getindex(itp::GriddedInterpolation{T,N,TCoefs,Interpolations.Gridded{Interpolations.Cubic{Interpolations.Natural}},K,P}, x::AbstractVector) where{T,N,TCoefs,K,P} + a,b,c,d = coefficients(itp.knots[1], itp.coefs) interpolate.(x, [a], [b], [c], [d], [itp.knots[1]]) end @@ -27,61 +27,97 @@ function interpolate(x, a, b, c, d, X, v=false) end +## TODO: Check which one of the below is faster. +function expand_array!(x2t, x2) + x2t[1] = x2[1] + x2t[end] = x2[end] + x2t[2:end-1] = cat(2, x2[2:end-1], x2[2:end-1])'[:] + x2t +end + +function expand_array2!(rhs, y) + n = length(y) + rhs[1] = y[1] # f₀(1) = y[1] + for i in 2:(n-1) + k = 2(i-1) # fᵢ = y[i] + rhs[k] = y[i] # fᵢ₊₁ = y[i] + rhs[k+1] = y[i] + end + rhs[n] = y[n] # fₙ = y[n] +end + +function expand_array(x2::AbstractArray{T,1}) where T + x2t = zeros(T,2(size(x2,1)-1)) + expand_array!(x2t, x2) +end + function coefficients(x, y; + force_linear_extrapolation = true, + boundary_condition = :natural) + A,rhs = spline_coef_equations(x, y, force_linear_extrapolation=force_linear_extrapolation, + boundary_condition=boundary_condition) + coefficients = A\rhs + a,b,c,d = [coefficients[n:4:end] for n in 1:4] + return a,b,c,d +end + +function spline_coef_equations(x, y; force_linear_extrapolation = true, boundary_condition = :natural ) - const x1 = 0. - const xend = 0. - - n = length(x) - A = zeros(n,n) - rx = zeros(n) - for i in 2:(n-1) - A[i-1,i] = 1.0/3.0*(x[i]-x[i-1]); - A[i,i] = 2.0/3.0*(x[i+1]-x[i-1]); - A[i+1,i] = 1.0/3.0*(x[i+1]-x[i]); - rx[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]) - end - - if boundary_condition == :natural - A[1,1] = 2.0 - A[2,1] = 0.0 - rx[1] = x1 - - A[n,n] = 2.0 - A[n-1,n] = 0.0 - rx[n] = xend - elseif boundary_condition == :periodic - A[1,1] = 2.0*(x[2]-x[1]) - A[2,1] = 1.0*(x[2]-x[1]) - rx[1] = 3.0*((y[2]-y[1])/(x[2]-x[1])-x1) - - A[n,n] = 2.0*(x[n]-x[n-1]) - A[n-1,n] =1.0*(x[n]-x[n-1]) - rx[n]=3.0*(xend-(y[n]-y[n-1])/(x[n]-x[n-1])) - else - error("Boundary condition $boundary_condition not recognised.") + valid_boundary_conditions = [:natural, :periodic, :notaknot, :quadratic] + if !(boundary_condition in valid_boundary_conditions) + error("Boundary Condition must be one of $valid_boundary_conditions not $boundary_condition") end - b = A \ rx; + n = length(x) + A = zeros(4(n-1),4(n-1)) + rhs = zeros(4(n-1)) + + # First 2(n-1) polynomials are of the form: + # fᵢ = aᵢx³ + bᵢx² +cᵢx + dᵢ + rhs[1:2(n-1)] = expand_array(y) - a = zeros(n) - c = zeros(n) - for i in 1:(n-1) - a[i]=1.0/3.0*(b[i+1]-b[i])/(x[i+1]-x[i]) - c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])- 1.0/3.0*(2.0*b[i]+b[i+1])*(x[i+1]-x[i]) + x2 = expand_array(x) + x2inds = expand_array(1:length(x)) + niind = 0 + for xi in 1:(length(x2)) + for ni in 1:4 + A[xi, (niind+ni)] = x2[xi]^(4-ni) + end + iseven(xi) && (niind += 4 ) end - if force_linear_extrapolation - b[1] = 0.0 - b[n] = 0.0 - a[n] = 0.0 - @show a[1] - h = x[n]-x[n-1]; - c[n] = 3a[n-1] * (h * h) + 2b[n-1] * h + c[n-1] + # Next polynomials from first derivative + # 3ax² + 2bx + c + 0 + niind = 1 + for xi in 2:(n-1) + rind = xi + 2(n-1) - 1 + A[rind,niind:niind+3] .= A[rind,(niind+4):(niind+7)] .= [(3*(x[xi]^2)),2x[xi],1,0] + A[rind,niind+4:niind+7] *= -1. + niind += 4 + end + + # Next polynomials from 2nd derivative + # 6ax + 2b + 0 + 0 + niind = 1 + for xi in 2:(n-1) + rind = xi + 2(n-1) + n - 3 + A[rind,niind:niind+3] .= A[rind,(niind+4):(niind+7)] .= [6*(x[xi]),2,0,0] + A[rind,niind+4:niind+7] *= -1. + niind += 4 end - return a, b, c + + # Next boundary conditions: + # Natural Spline Satisfies: + # 6a₁x₁ + 2b₁ = 0 + A[end-1,1] = 6x[1] + A[end-1,2] = 2. + # 6aₙxₙ₊₁ + 2bₙ = 0 + A[end,4(n-2)+1] = 6x[end] + A[end,4(n-2)+2] = 2. + + return A,rhs end #####