[WIP] A generalized PeriodicEuclidean (allows specifying the periods as a parallelogon)

.gitignore

@@ -2,5 +2,6 @@
 *.jl.*.cov
 *.jl.mem
 *.ji
+.vscode/
 benchmark/params.json
 Manifest.toml

src/metrics.jl

@@ -394,6 +394,20 @@ weuclidean(a, b, w) = WeightedEuclidean(w)(a, b)
     s3 = min(s2, p - s2)
     abs2(s3)
 end
+
+function (pe::PeriodicEuclidean{<:AbstractMatrix})(a::AbstractVector, b::AbstractVector)
+    p = parameters(pe)
+    ndim = length(a)
+    @assert size(p,1) == size(p, 2) == length(b) == ndim  # check dimension
+    s = p \ (a - b)      # decompose to p-basis, a - b = sx x⃗ + sy y⃗ + ...
+    mindistance = Inf
+    s = mod.(s, 1)       # move to first "quadrant"
+    @inbounds for k = 0:1<<ndim-1  # 2^{# of dimensions} possible choices
+        loc = sum(i->view(p,:,i) .* (k>>(i-1) & 1 == 0 ? s[i] : s[i]-1), 1:ndim)
+        mindistance = min(norm(loc), mindistance)
+    end
+    return mindistance
+end
 eval_end(::PeriodicEuclidean, s) = sqrt(s)
 peuclidean(a, b, p) = PeriodicEuclidean(p)(a, b)
 

test/test_dists.jl

@@ -982,3 +982,11 @@ end
     end
 end
 =#
+
+@testset "periodic - non-rectangular" begin
+    # basis vectors: [10, 10], [0, 10]
+    p = PeriodicEuclidean([10 0.0; 10.0 10.0])  # each column is a basis vector.
+    @test p([10.1, 10.2], [-10.0, 0.0]) ≈ sqrt(0.1^2 + 0.2^2)
+    @test p([9.9, 10.2], [-10.0, 0.0]) ≈ sqrt(0.1^2 + 0.2^2)
+    @test p([9.9, 9.8], [-10.0, 0.0]) ≈ sqrt(0.1^2 + 0.2^2)
+end