From f931a8ac5b3996aeb49f935403e68d33ca2bdb4a Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Fri, 8 Dec 2023 12:22:44 -0800 Subject: [PATCH 1/2] Optimization changes Some simple: - Remove unnecessary array allocation - Improve performance of _G function Some complicated - Avoid returning multiple function calls - Use let function to avoid captured variables (see https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured) --- src/GGDUtils.jl | 2 ++ src/interpolations.jl | 63 ++++++++++++++++++++++--------------------- src/subset_tools.jl | 20 +++++++------- 3 files changed, 44 insertions(+), 41 deletions(-) diff --git a/src/GGDUtils.jl b/src/GGDUtils.jl index 67c639c..6c7fb25 100644 --- a/src/GGDUtils.jl +++ b/src/GGDUtils.jl @@ -2,6 +2,8 @@ module GGDUtils using OMAS: OMAS +const inv_16pi = 1.0 / (16π) + export project_prop_on_subset! export get_subset_centers export get_prop_with_grid_subset_index diff --git a/src/interpolations.jl b/src/interpolations.jl index f8992ef..ca1be38 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -68,11 +68,9 @@ minimizing bending energy of a surface. http://www.geometrictools.com/Documentation/ThinPlateSplines.pdf Eq(28) """ function _G(x1::Tuple{U, U}, x2::Tuple{U, U}) where {U <: Real} - r = sqrt(sum((x1 .- x2) .^ 2)) - if r == 0 - return 0 - end - return r^2 * log(r) / 8 / π + r2 = sum((x1 .- x2) .^ 2) # r² + # Note this uses log(r) / 8π = log(r²) / 16π + return (r2 == 0) ? zero(U) : r2 * log(r2) * inv_16pi end """ @@ -87,21 +85,16 @@ minimum absolute value. This is done to avoid numerical issues with interpolatio Return values are conditioned y and inverse conditioning function. """ function _condition_y(y::Vector{T}) where {T <: Real} - do_log = false ylims = extrema(y) - if prod(ylims) > 0 - if ylims[2] / ylims[1] > 100 - do_log = true - end - end - if do_log - norm_by = minimum(abs.(ylims)) * sign(ylims[1]) - return log10.(y ./ norm_by), (cy) -> (10 .^ (cy)) * norm_by - else - norm_by = ylims[2] - ylims[1] - mean_y = mean(y) - return (y .- mean_y) ./ norm_by, (cy) -> (cy .* norm_by) .+ mean_y + do_log = (prod(ylims) > 0) && (ylims[2] / ylims[1] > 100) + norm_by = do_log ? minimum(abs.(ylims)) * sign(ylims[1]) : ylims[2] - ylims[1] + mean_y = mean(y) + cy = do_log ? log10.(y ./ norm_by) : (y .- mean_y) ./ norm_by + inv_cy = let norm_by = norm_by, do_log = do_log, mean_y = mean_y + x -> do_log ? (10.0 ^ x) * norm_by : (x * norm_by) + mean_y end + return cy, inv_cy + end """ @@ -129,6 +122,12 @@ function get_TPS_mats(x::Vector{Tuple{U, U}}) where {U <: Real} return Minv, N, (N' * Minv * N)^(-1) * N' * Minv, x end +function get_interp_val(r, z, x, a, b, inv_cy) + tot = sum(a[k] * _G((r, z), x[k]) for k in eachindex(a)) + tot += b[1] + r * b[2] + z * b[3] + return inv_cy(tot) +end + """ interp( y::Vector{T}, @@ -147,11 +146,10 @@ function interp( # From Eq(31) b = y2b * cy a = Minv * (cy - N * b) - function get_interp_val(r::Real, z::Real) - return inv_cy(sum(a .* [_G((r, z), xi) for xi ∈ x]) + sum(b .* [1, r, z])) + g = let x = x, a = a, b = b, inv_cy = inv_cy + (r, z) -> get_interp_val(r, z, x, a, b, inv_cy) end - get_interp_val(gp::Tuple{V, V}) where {V <: Real} = get_interp_val(gp...) - return get_interp_val + return g end """ @@ -322,10 +320,11 @@ function interp( prop_arr::Vector{T}, TPS_mats::Tuple{Matrix{U}, Matrix{U}, Matrix{U}, Vector{Tuple{U, U}}}, grid_subset_index::Int, - value_field::Symbol=:values, -) where {T <: edge_profiles__prop_on_subset, U <: Real} + value_field::Val{V}=Val(:values), +) where {T <: edge_profiles__prop_on_subset, U <: Real, V} prop = get_prop_with_grid_subset_index(prop_arr, grid_subset_index) - return interp(getfield(prop, value_field), TPS_mats) + field = getfield(prop, V) + return interp(field, TPS_mats) end const RHO_EXT_POS = [1.0001, 1.1, 5] @@ -363,9 +362,10 @@ function interp(eqt::OMAS.equilibrium__time_slice) prepend!(rhon_eq_ext, RHO_EXT_NEG) rz2psin = linear_interpolation((r_eq, z_eq), psinrz) psin2rhon = linear_interpolation(psin_eq_ext, rhon_eq_ext) - get_interp_val(r::Real, z::Real) = psin2rhon(rz2psin(r, z)) - get_interp_val(rz::Tuple{Real, Real}) = get_interp_val(rz...) - return get_interp_val + g = let psin2rhon=psin2rhon, rz2psin=rz2psin + (r, z) -> psin2rhon(rz2psin(r, z)) + end + return g end """ @@ -423,9 +423,10 @@ function interp( rz2rho::Function, ) where {T <: Real} itp = interp(prop, prof) - get_interp_val(r::Real, z::Real) = itp.(rz2rho(r, z)) - get_interp_val(rz::Tuple{Real, Real}) = get_interp_val(rz...) - return get_interp_val + g = let itp=itp + (r, z) -> itp(rz2rho(r, z)) + end + return g end """ diff --git a/src/subset_tools.jl b/src/subset_tools.jl index 2ff175f..1638574 100644 --- a/src/subset_tools.jl +++ b/src/subset_tools.jl @@ -96,6 +96,10 @@ function get_grid_subset_with_index( return subset end end + # BCL 12/8: Creates type instability, but maybe okay since it's a "simple" Union + # subset::Union{Int64, OMAS.edge_profiles__grid_ggd___grid_subset} + # + # Better would be to immediately throw an error or return nothing return 0 # Indicates failure end @@ -202,7 +206,7 @@ function get_subset_centers(space::OMAS.edge_profiles__grid_ggd___space, subset_space = get_subset_space(space, subset.element) grid_nodes = space.objects_per_dimension[1].object return [ - Tuple(mean([grid_nodes[node].geometry for node ∈ obj.nodes])) for + Tuple(mean(SVector{2}(grid_nodes[node].geometry) for node ∈ obj.nodes)) for obj ∈ subset_space ] end @@ -213,7 +217,7 @@ end from_subset::OMAS.edge_profiles__grid_ggd___grid_subset, to_subset::OMAS.edge_profiles__grid_ggd___grid_subset, space::OMAS.edge_profiles__grid_ggd___space, -) +) This function can be used to add another instance on a property vector representing the value in a new subset that can be taken as a projection from an existing larger subset. @@ -389,21 +393,17 @@ function Base.:∈( count = 0 for ele ∈ subset_bnd.element edge = edges[ele.object[1].index] - r_max = maximum([nodes[node].geometry[1] for node ∈ edge.nodes]) - r_min = minimum([nodes[node].geometry[1] for node ∈ edge.nodes]) + r_max = maximum(nodes[node].geometry[1] for node ∈ edge.nodes) + r_min = minimum(nodes[node].geometry[1] for node ∈ edge.nodes) if r_min <= r < r_max - z_max = maximum([nodes[node].geometry[2] for node ∈ edge.nodes]) + z_max = maximum(nodes[node].geometry[2] for node ∈ edge.nodes) if z < z_max count += 1 end end end # If it is even, the point is outside the boundary - if count % 2 == 1 - return true - else - return false - end + return count % 2 == 1 end function get_prop_with_grid_subset_index( From 4e4e9ba0acd8b456a82f6e0dc15451ac6b041151 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Thu, 21 Dec 2023 10:00:03 -0800 Subject: [PATCH 2/2] Lowest level interp should return overloaded func for (r,z) & ((r,z)) --- src/interpolations.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/interpolations.jl b/src/interpolations.jl index ca1be38..1d4f9c2 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -87,14 +87,13 @@ Return values are conditioned y and inverse conditioning function. function _condition_y(y::Vector{T}) where {T <: Real} ylims = extrema(y) do_log = (prod(ylims) > 0) && (ylims[2] / ylims[1] > 100) - norm_by = do_log ? minimum(abs.(ylims)) * sign(ylims[1]) : ylims[2] - ylims[1] + norm_by = do_log ? minimum(abs.(ylims)) * sign(ylims[1]) : ylims[2] - ylims[1] mean_y = mean(y) cy = do_log ? log10.(y ./ norm_by) : (y .- mean_y) ./ norm_by inv_cy = let norm_by = norm_by, do_log = do_log, mean_y = mean_y - x -> do_log ? (10.0 ^ x) * norm_by : (x * norm_by) + mean_y + x -> do_log ? (10.0^x) * norm_by : (x * norm_by) + mean_y end return cy, inv_cy - end """ @@ -123,7 +122,7 @@ function get_TPS_mats(x::Vector{Tuple{U, U}}) where {U <: Real} end function get_interp_val(r, z, x, a, b, inv_cy) - tot = sum(a[k] * _G((r, z), x[k]) for k in eachindex(a)) + tot = sum(a[k] * _G((r, z), x[k]) for k ∈ eachindex(a)) tot += b[1] + r * b[2] + z * b[3] return inv_cy(tot) end @@ -146,10 +145,10 @@ function interp( # From Eq(31) b = y2b * cy a = Minv * (cy - N * b) - g = let x = x, a = a, b = b, inv_cy = inv_cy - (r, z) -> get_interp_val(r, z, x, a, b, inv_cy) + return let x = x, a = a, b = b, inv_cy = inv_cy + g(r, z) = get_interp_val(r, z, x, a, b, inv_cy) + g(gp::Tuple{V, V}) where {V <: Real} = g(gp...) end - return g end """ @@ -362,7 +361,7 @@ function interp(eqt::OMAS.equilibrium__time_slice) prepend!(rhon_eq_ext, RHO_EXT_NEG) rz2psin = linear_interpolation((r_eq, z_eq), psinrz) psin2rhon = linear_interpolation(psin_eq_ext, rhon_eq_ext) - g = let psin2rhon=psin2rhon, rz2psin=rz2psin + g = let psin2rhon = psin2rhon, rz2psin = rz2psin (r, z) -> psin2rhon(rz2psin(r, z)) end return g @@ -423,7 +422,7 @@ function interp( rz2rho::Function, ) where {T <: Real} itp = interp(prop, prof) - g = let itp=itp + g = let itp = itp (r, z) -> itp(rz2rho(r, z)) end return g