From e34e9c3df20e68cb4aece99386e6c27718cb4dfa Mon Sep 17 00:00:00 2001 From: cristianeferreira Date: Mon, 26 Aug 2024 12:00:24 -0300 Subject: [PATCH] Remove matrix division from cubic splines (#19) * Remove matrix division * Include blank line at the end --- src/splines.jl | 141 ++++++++++++++++--------------------------------- 1 file changed, 45 insertions(+), 96 deletions(-) diff --git a/src/splines.jl b/src/splines.jl index 0a2d1bc..e6af44e 100644 --- a/src/splines.jl +++ b/src/splines.jl @@ -14,36 +14,34 @@ mutable struct Spline{T} end end -# Given a polinomial index, returns the indexer of the last parameter before the first parameter of the referenced polinomial -_base_param_index_(poly_index::Int) = (poly_index-1)*4 - # Performs natural cubic spline interpolation function splineint(s::Spline{T}, x_out::Number) where {T<:Number} - local poly_index::Int = 1 - local base_idx::Int - + if x_out > s.x[end] # Extrapolation after last point - poly_index = length(s.x) - 1 - base_idx = _base_param_index_(poly_index) - return (x_out - s.x[end])*(s.params[base_idx + 2] + 2*s.params[base_idx + 3]*s.x[end] + 3*s.params[base_idx + 4]*(s.x[end]^2)) + s.y[end] + n = length(s.x) + x = x_out - s.x[n] + dxn = s.x[n] - s.x[n-1] + b = s.params[n-3] + 2*s.params[n-2]*dxn + 3*s.params[n-1]*dxn*dxn + return s.y[n] + b*x + elseif x_out < s.x[1] # Extrapolation before first point - base_idx = 0 - return (x_out - s.x[1])*(s.params[base_idx + 2] + 2*s.params[base_idx + 3]*s.x[1] + 3*s.params[base_idx + 4]*(s.x[1]^2)) + s.y[1] + return s.params[1] + s.params[2]*(x_out - s.x[1]) + else - # Interplation - while x_out > s.x[poly_index+1] - poly_index += 1 + # Find polynomial + i = 1 + while x_out > s.x[i+1] + i = i + 1 end - - base_idx = _base_param_index_(poly_index) - - #P1 P2 P3 ... - #1 2 3 4 5 6 7 8 9 10 11 12 - #a + bx + cx^2 + dx^3 a + bx + cx^2 + dx^3 a + bx + cx^2 + dx^3 ... - return s.params[base_idx + 1] + s.params[base_idx + 2]*x_out + - s.params[base_idx + 3]*(x_out^2) + s.params[base_idx + 4]*(x_out^3) + + x = x_out - s.x[i] + a = s.params[4*i - 3] + b = s.params[4*i - 2] + c = s.params[4*i - 1] + d = s.params[4*i] + return a + b*x + c*x*x + d*x*x*x end end @@ -59,81 +57,32 @@ end # Build a Spline object by fitting 3rd order polynomials around points (x_in, y_in) function splinefit(x_in::Vector{T}, y_in::Vector{Float64}) :: Spline{T} where {T} - # - # TODO: optimize. See http://www.math.ntnu.no/emner/TMA4215/2008h/cubicsplines.pdf - # - points_count = length(x_in) - @assert points_count == length(y_in) "x_in and y_in doesn't conform on sizes." - matrix_n = 4*(points_count-1) # the main matrix is a square matrix matrix_n by matrix_n - - A = zeros(matrix_n, matrix_n) - b = zeros(matrix_n) - - # Known values for polys - - # First Point - A[1,1] = 1.0 - A[1,2] = x_in[1] - A[1,3] = x_in[1]^2 - A[1,4] = x_in[1]^3 - b[1] = y_in[1] - - # Last Point - base_idx = 4*(points_count - 2) - A[2, base_idx + 1] = 1.0 - A[2, base_idx + 2] = x_in[points_count] - A[2, base_idx + 3] = x_in[points_count]^2 - A[2, base_idx + 4] = x_in[points_count]^3 - b[2] = y_in[points_count] - - # Inner Points - row = 3 - for i in 2:(points_count-1) - # Connecting to left poly - A[row, (i-2)*4 + 1] = 1.0 - A[row, (i-2)*4 + 2] = x_in[i] - A[row, (i-2)*4 + 3] = x_in[i]^2 - A[row, (i-2)*4 + 4] = x_in[i]^3 - b[row] = y_in[i] - - row += 1 - - # Connecting to right poly - A[row, (i-1)*4 + 1] = 1.0 - A[row, (i-1)*4 + 2] = x_in[i] - A[row, (i-1)*4 + 3] = x_in[i]^2 - A[row, (i-1)*4 + 4] = x_in[i]^3 - b[row] = y_in[i] - - row += 1 - - # Conditions on first order derivatives - A[row, (i-2)*4 + 2] = 1.0 - A[row, (i-2)*4 + 3] = 2.0*x_in[i] - A[row, (i-2)*4 + 4] = 3.0*(x_in[i]^2) - A[row, (i-1)*4 + 2] = -1.0 - A[row, (i-1)*4 + 3] = -2.0*x_in[i] - A[row, (i-1)*4 + 4] = -3.0*(x_in[i]^2) - - row += 1 - - # Conditions on second order derivatives - A[row, (i-2)*4 + 3] = 2.0 - A[row, (i-2)*4 + 4] = 6.0 * x_in[i] - A[row, (i-1)*4 + 3] = -2.0 - A[row, (i-1)*4 + 4] = -6.0 * x_in[i] - - row += 1 + n = length(x_in) + @assert n == length(y_in) "x_in and y_in doesn't conform on sizes." + + H = [X[i+1] - X[i] for i in 1:(n-1)] # i in [0,...,n-2] + A = Y + Alpha = [(3/H[i])*(A[i+1] - A[i]) - (3/H[i-1])*(A[i] - A[i-1]) for i in 2:(n-1)] + Alpha = [0; Alpha] + + L = ones(n) + Mu = zeros(n) + Z = zeros(n) + B = zeros(n) + C = zeros(n) + D = zeros(n) + + for i in 2:n-1 + L[i] = 2 * (X[i+1] - X[i-1]) - H[i-1] * Mu[i-1] + Mu[i] = H[i] / L[i] + Z[i] = (Alpha[i] - H[i-1]*Z[i-1]) / L[i] end - # Conditions for natural cubic spline - A[row, 3] = 2.0 - A[row, 4] = 6.0 * x_in[1] - - row += 1 - - A[row, (points_count-2)*4 + 3] = 2.0 - A[row, (points_count-2)*4 + 4] = 6.0 * x_in[points_count] + for i in n-1:-1:1 + C[i] = Z[i] - Mu[i] * C[i+1] + B[i] = (A[i+1] - A[i]) / H[i] - H[i] * (C[i+1] + 2*C[i]) / 3 + D[i] = (C[i+1] - C[i]) / (3*H[i]) + end - return Spline(x_in, y_in, A \ b) + return reduce(vcat, [[A[i], B[i], C[i], D[i]] for i in 1:n-1]) end