add interpolate(kp; length, density) kwarg methods (cf. JuliaMolSim…

…/DFTK.jl#496)

- bump version to 0.4.1

Project.toml

@@ -1,7 +1,7 @@
 name = "Brillouin"
 uuid = "23470ee3-d0df-4052-8b1a-8cbd6363e7f0"
 authors = ["Thomas Christensen <[email protected]> and contributors"]
-version = "0.4.0"
+version = "0.4.1"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

src/interpolate-paths.jl

@@ -89,12 +89,18 @@ function latticize!(kpi::KPathInterpolant)
 end
 
 """
-    interpolate(kp::KPath, N::Integer) --> KPathInterpolant
+    interpolate(kp::KPath, length::Integer)
+    interpolate(kp::KPath; length::Integer, density::Real) --> KPathInterpolant
 
-Return an interpolant of `kp` with approximately `N` points distributed approximately
+Return an interpolant of `kp` with approximately `length` points distributed approximately
 equidistantly across the full **k**-path (equidistance is measured in a Cartesian metric).
 
-Note that the interpolant may contain fewer or more points than `N` (typically fewer). 
+Note that the interpolant may contain fewer or more points than `length` (typically fewer).
+`length` can also be provided as a keyword argument.
+
+As an alternative to specifying the desired total number of interpolate points via `length`,
+a desired density per unit (reciprocal) length can be specified via the keyword argument
+`density`.
 
 See also [`interpolate(::AbstractVector{::AbstractVector{<:Real}}, ::Integer)`](@ref).
 """
@@ -118,7 +124,54 @@ function interpolate(kp::KPath{D}, N::Integer) where D
                                              length=Nᵢ)[1:end-1]))
             push!(labels[j], length(kipaths[j])+1 => path[i+1])
         end
-        push!(kipaths[j], points(kp)[path[end]])
+        push!(kipaths[j], points(kp)[last(path)])
+    end
+
+    return KPathInterpolant(kipaths, labels, basis(kp), kp.basisenum)
+end
+
+function interpolate(kp::KPath;
+            length::Union{Integer, Nothing}=nothing,
+            density::Union{Real, Nothing}=nothing)
+    if length === nothing && density === nothing
+        throw(ArgumentError("must set either `length` or `density` keyword arguments explicitly or use positional argument for `length`"))
+    elseif length !== nothing && density !== nothing
+        throw(ArgumentError("cannot define `length` and `density` simultaneously"))
+    elseif length !== nothing && density === nothing
+        return interpolate(kp, length)
+    elseif length === nothing && density !== nothing
+        return interpolate_via_density(kp, density)
+    end
+
+    error("unreachable reached; please file an issue")
+end
+
+function interpolate_via_density(kp::KPath{D}, density::Real) where D
+    kpᶜ = kp.basisenum[] === CARTESIAN ? kp : cartesianize(kp)
+
+    kipaths = [Vector{SVector{D, Float64}}() for _ in 1:length(paths(kp))]
+    labels = [Dict{Int, Symbol}() for _ in 1:length(paths(kp))]
+    for (j, path) in enumerate(paths(kpᶜ)) # over connected paths
+        push!(labels[j], 1 => first(path))
+        for i in 1:length(path)-1 # over linear segments
+            klabᵢ   = path[i]
+            klabᵢ₊₁ = path[i+1]
+
+            # compute required number of points for current segment
+            kᶜᵢ   = points(kpᶜ)[klabᵢ]   # start
+            kᶜᵢ₊₁ = points(kpᶜ)[klabᵢ₊₁] # end
+            distᵢ = norm(kᶜᵢ₊₁ - kᶜᵢ)
+            Nᵢ    = max(ceil(Int, distᵢ*density), 2)
+
+            # compute interpolation of current segment (excluding end point)
+            kᵢ    = points(kp)[klabᵢ]
+            kᵢ₊₁  = points(kp)[klabᵢ₊₁]
+            ks    = range(kᵢ, kᵢ₊₁, length=Nᵢ)[1:end-1]
+
+            append!(kipaths[j], ks)
+            push!(labels[j], length(kipaths[j])+1 => klabᵢ₊₁)
+        end
+        push!(kipaths[j], points(kp)[last(path)]) # add previously omitted end point
     end
 
     return KPathInterpolant(kipaths, labels, basis(kp), kp.basisenum)
@@ -143,7 +196,7 @@ function splice(kp::KPath{D}, N::Integer) where D
                                              length=N+2)[1:end-1]))
             push!(labels[j], length(kipaths[j])+1 => path[i+1])
         end
-        push!(kipaths[j], points(kp)[path[end]])
+        push!(kipaths[j], points(kp)[last(path)])
     end
 
     return KPathInterpolant(kipaths, labels, kp.basis, kp.basisenum)

test/kpaths.jl

@@ -1,7 +1,7 @@
 import Crystalline
 using Brillouin.KPaths: cartesianize, latticize
 
-@testset "KPaths" begin
+@testset "KPath & KPathInterpolant" begin
     # --- `cumdist` ---
     kvs_2d = [[0,0], [0,1], [1,1], [1,0], [0,0], [1,1], [-1,-2], [2,-1]]
     @test cumdists(kvs_2d) ≈ [0,1,2,3,4,4+√2,4+√2+√13, 4+√2+√13+√10]
@@ -97,6 +97,7 @@ using Brillouin.KPaths: cartesianize, latticize
     Rs′ = Crystalline.DirectBasis([1, 0, 0], [0.3, 0.8, 0], [-1.6, 0.8, 0.9]) # neither all-obtuse nor all-acute
     @test_throws DomainError Brillouin.KPaths.extended_bravais(1, "aP", Rs′, Val(3))       # `_throw_basis_required`
 
+
     # --- `KPathInterpolant` ---
     # `interpolate`
     N = 100
@@ -105,6 +106,7 @@ using Brillouin.KPaths: cartesianize, latticize
 
     # all points in `kp` must be explicitly contained in `kpi`
     @test all(∈(kpi), values(points(kp)))
+    @test kpi == interpolate(kp, length=N)
 
     # `splice`
     kps = splice(kp, N)
@@ -126,4 +128,14 @@ using Brillouin.KPaths: cartesianize, latticize
     @test typeof(kpi) === typeof(kpi′) === KPathInterpolant{2}
     @test kpi′ ≈ kpi # test that `interpolate` commutes with `latticize`/`cartesianize`
     @test Brillouin.cartesianize!(latticize!(kpi)) ≈ kpi
-end
+
+    # interpolate via `density` kwarg
+    for ρ = [1, 10, 20, 21, 100, 3.0]
+        kpi_ρ = interpolate(kp, density=10)
+        @test kpi_ρ ≈ interpolate(kp, length=length(kpi_ρ))
+    end
+
+    # interpolate with wrong kwargs
+    @test_throws ArgumentError interpolate(kp)
+    @test_throws ArgumentError interpolate(kp; length=10, density=2)
+end # testset