Skip to content

Conversation

@hervasa2
Copy link
Collaborator

Continuing the discussion on #545 @fhagemann suggested to switch to FritschCarlsonMonotonicInterpolation.
This provides an elegant solution which is built into Interpolations.

Copy link
Collaborator

@fhagemann fhagemann left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good to me.

For reference, this is what the examples in the previous PR #545 look like now:

z = [1,5,6,10]
ρ = [1,2,3,4]
scatter(z, ρ)
plot!(range(extrema(z)..., length = 100),z ->  spline(z))
image
z = [1,5,6.5,10]
ρ = [1,2,3,4]
ρ = [1,2,3,2]
scatter(z, ρ)
plot!(range(extrema(z)..., length = 100),z ->  spline(z))
image
z = [1,5,7.6,10]
ρ = [1,2,3,2]
scatter(z, ρ)
plot!(range(extrema(z)..., length = 100),z ->  spline(z))
image

@fhagemann fhagemann added enhancement Improvement of existing features convenience Improve user-friendliness labels Oct 26, 2025
@fhagemann
Copy link
Collaborator

Also, seems like we can decrease the bound for Interpolations back to 0.14 then, as we don't require the use of cubic_spline_interpolation anymore.

@fhagemann
Copy link
Collaborator

The only thing that does not (yet) work with this PR: FritschCarlsonMonotonicInterpolation does not seem to extrapolate, so if the user were to evaluate this impurity density outside the z range where it is defined, this method will result in a wordy error message:

z = [1,5,6,10]
ρ = [1,2,3,4]
plot!(0:1e-3:11, z ->  spline(z))
ERROR: BoundsError: attempt to access 4-element interpolate(::SVector{4, Float64}, ::SVector{4, Float64}, FritschCarlsonMonotonicInterpolation()) with element type Float64 with indices SOneTo(4) at index [0.0]
Stacktrace:
  [1] throw_boundserror(A::Interpolations.MonotonicInterpolation{…}, I::Tuple{…})
    @ Base ./essentials.jl:14
  [2] (::Interpolations.MonotonicInterpolation{…})(x::Float64)
    @ Interpolations ~/.julia/packages/Interpolations/dR5oF/src/monotonic/monotonic.jl:211
  [3] (::var"#1#2")(z::Float64)
    @ Main ./REPL[16]:1
  [4] iterate
    @ ./generator.jl:48 [inlined]
  [5] _collect(c::StepRangeLen{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
    @ Base ./array.jl:811
  [6] collect_similar(cont::StepRangeLen{…}, itr::Base.Generator{…})
    @ Base ./array.jl:720
  [7] map(f::Function, A::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64})
    @ Base ./abstractarray.jl:3371
  [8] _compute_y
    @ ~/.julia/packages/RecipesPipeline/BGM3l/src/series.jl:72 [inlined]
  [9] _compute_xyz(x::StepRangeLen{…}, y::Function, z::Nothing, nice_error::Bool)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/series.jl:86
 [10] macro expansion
    @ ~/.julia/packages/RecipesPipeline/BGM3l/src/series.jl:140 [inlined]
 [11] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, ::Type{RecipesPipeline.SliceIt}, x::Any, y::Any, z::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesBase/BRe07/src/RecipesBase.jl:300
 [12] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/user_recipe.jl:38
 [13] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/RecipesPipeline.jl:72
 [14] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
    @ Plots ~/.julia/packages/Plots/8ZnR3/src/plot.jl:223
 [15] plot!(::Plots.Plot, ::Any, ::Vararg{Any}; kw...)
    @ Plots ~/.julia/packages/Plots/8ZnR3/src/plot.jl:213
 [16] plot!(::Plots.Plot, ::Any, ::Any)
    @ Plots ~/.julia/packages/Plots/8ZnR3/src/plot.jl:208
 [17] plot!(::Any, ::Vararg{Any}; kw...)
    @ Plots ~/.julia/packages/Plots/8ZnR3/src/plot.jl:202
 [18] plot!(::Any, ::Any)
    @ Plots ~/.julia/packages/Plots/8ZnR3/src/plot.jl:194
 [19] top-level scope
    @ REPL[16]:1
Some type information was truncated. Use `show(err)` to see complete types.

Should we add some catch for this (e.g. if z < zmin, just take the first value of ρ and if z > zmax, take the last value of ρ), throw an explicit warning or just leave this as is?

@hervasa2
Copy link
Collaborator Author

Should we add some catch for this (e.g. if z < zmin, just take the first value of ρ and if z > zmax, take the last value of ρ), throw an explicit warning or just leave this as is?

With the cubic spline it was the same story. But I would put the check in get_impurity_density. If the user use idm.spline manually there should be no check.

Copy link
Collaborator

@fhagemann fhagemann left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perfect, thanks!

@fhagemann
Copy link
Collaborator

Tests pass locally, I am going to merge this now.

@fhagemann fhagemann merged commit 7336b54 into JuliaPhysics:main Oct 26, 2025
9 checks passed
fhagemann added a commit that referenced this pull request Oct 26, 2025
…olation` (#547)

* switch to FritschCarlsonMonotonicInterpolation

* check that point is in boule

* assert for SplineBouleImpurityDensity

* Lower compat for Interpolations back to 0.14

* Update boule impurity tests

---------

Co-authored-by: Felix Hagemann <[email protected]>
(cherry picked from commit 7336b54)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

convenience Improve user-friendliness enhancement Improvement of existing features

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants