From c72bfd80d315837d4303ed7d468dd4057ddbcaef Mon Sep 17 00:00:00 2001 From: RJDennis Date: Mon, 24 May 2021 18:05:28 +0100 Subject: [PATCH] Replaced generated function --- .travis.yml | 8 +-- Project.toml | 2 +- src/piecewise_linear.jl | 154 +++++++++++++++++++--------------------- test/runtests.jl | 2 +- 4 files changed, 79 insertions(+), 87 deletions(-) diff --git a/.travis.yml b/.travis.yml index adc415e..87e9e25 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,12 +1,12 @@ language: julia julia: - - 1.1 - - 1.2 - - 1.3 + - 1.5 + - 1.6 + - nightly #script: # - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi -# - julia -e 'Pkg.clone(pwd()); Pkg.test("SmolyakApprox", coverage=true)' +# - julia -e 'Pkg.clone(pwd()); Pkg.test("PiecewiseLinearApprox", coverage=true)' after_success: # push coverage results to Codecov - julia -e 'cd(Pkg.dir("PiecewiseLinearApprox")); Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())' diff --git a/Project.toml b/Project.toml index 8b4e89e..0f4834c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PiecewiseLinearApprox" uuid = "fcfa6960-8e2e-11e9-2cde-29dd09221fb3" authors = ["Richard Dennis "] -version = "0.1.1" +version = "0.1.2" [compat] julia = "^1.1" diff --git a/src/piecewise_linear.jl b/src/piecewise_linear.jl index 001fd48..4a5c054 100644 --- a/src/piecewise_linear.jl +++ b/src/piecewise_linear.jl @@ -7,8 +7,29 @@ of dimensions function piecewise_linear_nodes(n::S,domain = [1.0,-1.0]) where {S <: Integer} - nodes = collect(range(domain[2],domain[1],length=n)) # The nodes are ordered from lowest to highest + if n == 1 + + return [(domain[1]+domain[2])/2] + + elseif n == 2 + + return [domain[2], domain[1]] + + else + + nodes = zeros(n) + nodes[1] = domain[2] + nodes[n] = domain[1] + + inc = (domain[1]-domain[2])/(n-1) + for i = 2:div(n,2) + nodes[i] = domain[2] + (i-1)*inc + nodes[n-i+1] = domain[1] - (i-1)*inc + end + return nodes + + end end @@ -16,23 +37,33 @@ const linear_nodes = piecewise_linear_nodes function bracket_nodes(x::Array{T,1},point::T) where {T <: AbstractFloat} - if point <= x[1] - return (1,2) - elseif point >= x[end] - return (length(x) - 1, length(x)) - else - return (sum(point .> x), sum(point .> x) + 1) + if point <= x[1] + return (1,2) + elseif point >= x[end] + return (length(x) - 1, length(x)) + else + y = 0 + for i in x + if i < point + y += 1 + else + break + end end + return (y, y+1) + end end - + function bracket_nodes(x::Union{NTuple{N,Array{T,1}},Array{Array{T,1},1}},point::Array{T,1}) where {T <: AbstractFloat, N} - bracketing_nodes = Array{Int64}(undef,2,length(x)) - for i = 1:length(point) - bracketing_nodes[:,i] .= bracket_nodes(x[i],point[i]) - end - return bracketing_nodes + bracketing_nodes = Array{Int64,2}(undef,2,length(point)) + for i = 1:length(point) + bracketing_nodes[:,i] .= bracket_nodes(x[i],point[i]) + end + + return bracketing_nodes + end function piecewise_linear_weight(x::Array{T,1},point::T) where {T <: AbstractFloat} @@ -46,7 +77,7 @@ end function piecewise_linear_weights(x::Union{NTuple{N,Array{T,1}},Array{Array{T,1},1}},point::Array{T,1}) where {T <: AbstractFloat, N} - weights = Array{T}(undef,length(point)) + weights = Array{T,1}(undef,length(point)) for i = 1:length(point) weights[i] = piecewise_linear_weight(x[i],point[i]) end @@ -72,12 +103,12 @@ end function select_relevant_data(y::AbstractArray{T,N},bracketing_grid_points::Array{S,2}) where {T <: AbstractFloat, S <: Integer, N} - data = zeros(size(bracketing_grid_points,1)) - for i = 1:length(data) - data[i] = y[CartesianIndex(tuple(bracketing_grid_points[i,:]...))] - end + data = zeros(size(bracketing_grid_points,1)) + for i = 1:length(data) + data[i] = y[CartesianIndex(tuple(bracketing_grid_points[i,:]...))] + end - return data + return data end @@ -93,7 +124,7 @@ function piecewise_linear_evaluate(y::AbstractArray{T,N},x::Union{NTuple{N,Array for j = d:-1:1 - new_data = zeros(Int(length(data)/2)) + new_data = zeros(div(length(data),2)) for i = 1:length(new_data) new_data[i] = data[2*(i-1)+1] + w[j]*(data[2*i]-data[2*(i-1)+1]) end @@ -146,42 +177,7 @@ function piecewise_linear_evaluate(y::AbstractArray{T,N},x::Union{NTuple{N,Array end -#= - -function locate_point_below(x::Array{T,1},point::T) where {T <: AbstractFloat} - - if point <= x[1] - return 1 - elseif point >= x[end] - return length(x) - 1 - else - return sum(point .> x) - end - -end - -function piecewise_linear_evaluate(y::Array{T,N},x::Union{NTuple{N,Array{T,1}},Array{Array{T,1},1}},point::Union{T,Array{T,1}}) where {T <: AbstractFloat, N} - - if length(x) == 1 - x_below = locate_point_below(x[1],point[1]) - x_above = x_below + 1 - y_hat = y[x_below] + ((point[1] - x[1][x_below])/(x[1][x_above] - x[1][x_below]))*(y[x_above] - y[x_below]) - return y_hat - else - x_below = locate_point_below(x[end],point[end]) - x_above = x_below + 1 - yy = zeros(size(y)[1:end-1]) - for i in CartesianIndices(yy) - yy[i] = y[CartesianIndex(i,x_below)] + ((point[end] - x[end][x_below])/(x[end][x_above] - x[end][x_below]))*(y[CartesianIndex(i,x_above)] - y[CartesianIndex(i,x_below)]) - end - x = x[1:end-1] - point = point[1:end-1] - piecewise_linear_evaluate(yy,x,point) - end - -end - -=# +# Functions to handle the 1-D case function piecewise_linear_evaluate(y::Array{T,1},x::Array{T,1},point::T) where {T <: AbstractFloat} @@ -206,36 +202,32 @@ function piecewise_linear_evaluate(y::Array{T,1},nodes::Array{T,1}) where {T <: end -@generated function grid_reshape(f::Array{T,N},grid::NTuple{N,Array{T,1}}) where {T,N} +# Function to transform from one uniform grid to another - i_vars = Array{Symbol}(undef,N) - grid_vars = Array{String}(undef,N) - for i = 1:N - i_vars[i] = Symbol("i$i") - grid_vars[i] = "grid[$i][i$i]" - end +function grid_reshape(f::Array{T,N},grid::NTuple{N,Array{T,1}}) where {T<:AbstractFloat,N} - inner = :(y[$(i_vars...)] = piecewise_linear_evaluate(f,old_grid,[$(Meta.parse.(grid_vars)...)])) + #1. Construct the old grid + + old_grid = Array{Array{T,1},1}(undef,N) + for i = 1:N + piecewise_linear_nodes + old_grid[i] = piecewise_linear_nodes(size(f,i),[grid[i][end],grid[i][1]]) + end - outer = inner + #2. Initialize the new y - for i = N:-1:1 - outer = :( - for $(i_vars[i]) = 1:length(grid[$i]) - $outer - end - ) - end + y = Array{T,N}(undef,tuple(length.(grid)...)) - final = :(old_grid = Array{Array{T,1},1}(undef,N); - for i = 1:N; - old_grid[i] = [range(grid[i][1],grid[i][end],length=size(f,i));]; - end; - y = Array{T,N}(undef,length.(grid)); - $outer; - return y - ) + #3. Interpolate to fill y + + point = Array{T,1}(undef,N) + for i in CartesianIndices(y) + for j = 1:N + point[j] = grid[j][i[j]] + end + y[i] = piecewise_linear_evaluate(f,old_grid,point) + end - return final + return y end diff --git a/test/runtests.jl b/test/runtests.jl index 46a524c..220274e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,7 +10,7 @@ using Test end - x1 = piecewise_linear_nodes(11,[3.0,1.0]) + x1 = piecewise_linear_nodes(11,[3.0,1.5]) x2 = piecewise_linear_nodes(21,[2.5,0.5]) x = (x1,x2)