From 27089556acbdcd4f19ef707dde7549b3b962bd82 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <joshua.lampert@uni-hamburg.de> Date: Wed, 18 Dec 2024 17:55:25 +0100 Subject: [PATCH 01/13] add parametric function for Chain --- src/geometries/polytopes.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/geometries/polytopes.jl b/src/geometries/polytopes.jl index ac16288f9..571306cd3 100644 --- a/src/geometries/polytopes.jl +++ b/src/geometries/polytopes.jl @@ -8,7 +8,7 @@ We say that a geometry is a K-polytope when it is a collection of "flat" sides that constitute a `K`-dimensional subspace. They are called chain, polygon and polyhedron respectively for 1D (`K=1`), 2D (`K=2`) and 3D (`K=3`) subspaces. -The parameter `K` is also known as the rank or parametric dimension +The parameter `K` is also known as the rank or parametric dimension of the polytope (<https://en.wikipedia.org/wiki/Abstract_polytope>). The term polytope expresses a particular combinatorial structure. A polyhedron, @@ -156,6 +156,16 @@ function angles(c::Chain) map(i -> ∠(vs[i - 1], vs[i], vs[i + 1]), i1:i2) end +function (c::Chain)(t) + if t < 0 || t > 1 + throw(DomainError(t, "c(t) is not defined for t outside [0, 1].")) + end + s = collect(segments(c)) + N = length(s) + k = max(1, ceil(Int, N * t)) + s[k](N * t - k + 1) +end + # implementations of Chain include("polytopes/segment.jl") include("polytopes/rope.jl") From 9c06204f385975711dfbe14f67815cd0b034952f Mon Sep 17 00:00:00 2001 From: Joshua Lampert <joshua.lampert@uni-hamburg.de> Date: Wed, 18 Dec 2024 18:15:48 +0100 Subject: [PATCH 02/13] add tests --- test/polytopes.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/polytopes.jl b/test/polytopes.jl index 2a50b8baa..137b72187 100644 --- a/test/polytopes.jl +++ b/test/polytopes.jl @@ -264,6 +264,22 @@ end r = Ring(merc.([(0, 0), (1, 0), (1, 1), (0, 1)])) @test crs(centroid(r)) === crs(r) + # parameterization + ri = Ring(latlon(45, 0), latlon(45, 90), latlon(90, 90)) + @test ri(T(0)) == latlon(45, 0) + @test ri(T(0.25)) == latlon(45, 67.5) + @test ri(T(0.5)) == latlon(67.5, 90) + @test ri(T(0.75)) == latlon(78.75, 67.5) + @test ri(T(1)) == latlon(45, 0) + + ro = Rope(latlon(45, 0), latlon(45, 90), latlon(90, 90)) + @test ro(T(0)) == latlon(45, 0) + @test ro(T(0.25)) == latlon(45, 45) + @test ro(T(0.5)) == latlon(45, 90) + @test ro(T(0.75)) == latlon(67.5, 90) + @test ro(T(1)) == latlon(90, 90) + + ri = Ring(cart.([(1, 1), (2, 2), (3, 3)])) ro = Rope(cart.([(1, 1), (2, 2), (3, 3)])) @test sprint(show, ri) == "Ring((x: 1.0 m, y: 1.0 m), (x: 2.0 m, y: 2.0 m), (x: 3.0 m, y: 3.0 m))" From 192b40ff830c784c7a11fdeae3e3fc1f20a64dd3 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Wed, 18 Dec 2024 18:18:16 +0100 Subject: [PATCH 03/13] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/polytopes.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/polytopes.jl b/test/polytopes.jl index 137b72187..9a66e7e1e 100644 --- a/test/polytopes.jl +++ b/test/polytopes.jl @@ -279,7 +279,6 @@ end @test ro(T(0.75)) == latlon(67.5, 90) @test ro(T(1)) == latlon(90, 90) - ri = Ring(cart.([(1, 1), (2, 2), (3, 3)])) ro = Rope(cart.([(1, 1), (2, 2), (3, 3)])) @test sprint(show, ri) == "Ring((x: 1.0 m, y: 1.0 m), (x: 2.0 m, y: 2.0 m), (x: 3.0 m, y: 3.0 m))" From 05a195be68055a2da3cf8bcd49a970aac5bc1017 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <joshua.lampert@uni-hamburg.de> Date: Wed, 18 Dec 2024 21:53:15 +0100 Subject: [PATCH 04/13] reduce allocations --- src/geometries/polytopes.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/geometries/polytopes.jl b/src/geometries/polytopes.jl index 571306cd3..7f642b684 100644 --- a/src/geometries/polytopes.jl +++ b/src/geometries/polytopes.jl @@ -160,10 +160,11 @@ function (c::Chain)(t) if t < 0 || t > 1 throw(DomainError(t, "c(t) is not defined for t outside [0, 1].")) end - s = collect(segments(c)) - N = length(s) - k = max(1, ceil(Int, N * t)) - s[k](N * t - k + 1) + v = vertices(c) + n = length(v) - !isclosed(c) + k = max(1, ceil(Int, n * t)) + s, _ = iterate(segments(c), k) + s(n * t - k + 1) end # implementations of Chain From dd528cc357500bdcd9a6bc428411062e906b0e69 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Wed, 18 Dec 2024 22:09:14 +0100 Subject: [PATCH 05/13] Update src/geometries/polytopes.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Júlio Hoffimann <julio.hoffimann@gmail.com> --- src/geometries/polytopes.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/geometries/polytopes.jl b/src/geometries/polytopes.jl index 7f642b684..80e2b0f98 100644 --- a/src/geometries/polytopes.jl +++ b/src/geometries/polytopes.jl @@ -160,8 +160,7 @@ function (c::Chain)(t) if t < 0 || t > 1 throw(DomainError(t, "c(t) is not defined for t outside [0, 1].")) end - v = vertices(c) - n = length(v) - !isclosed(c) + n = nvertices(c) - !isclosed(c) k = max(1, ceil(Int, n * t)) s, _ = iterate(segments(c), k) s(n * t - k + 1) From fbaaed58e4b8e582f666c916599f2f30ffcd3cdf Mon Sep 17 00:00:00 2001 From: Joshua Lampert <joshua.lampert@uni-hamburg.de> Date: Wed, 18 Dec 2024 22:35:26 +0100 Subject: [PATCH 06/13] fix --- src/geometries/polytopes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/geometries/polytopes.jl b/src/geometries/polytopes.jl index 80e2b0f98..c57cd7c17 100644 --- a/src/geometries/polytopes.jl +++ b/src/geometries/polytopes.jl @@ -161,9 +161,9 @@ function (c::Chain)(t) throw(DomainError(t, "c(t) is not defined for t outside [0, 1].")) end n = nvertices(c) - !isclosed(c) - k = max(1, ceil(Int, n * t)) + k = min(n - 1, floor(Int, n * t)) s, _ = iterate(segments(c), k) - s(n * t - k + 1) + s(n * t - k) end # implementations of Chain From 28d07dcd6da823ce7504952aeea023f314a357f6 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <joshua.lampert@uni-hamburg.de> Date: Thu, 19 Dec 2024 11:59:06 +0100 Subject: [PATCH 07/13] use measure --- src/geometries/polytopes.jl | 15 +++++++++++---- test/polytopes.jl | 12 ++++++------ 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/src/geometries/polytopes.jl b/src/geometries/polytopes.jl index c57cd7c17..bd922408e 100644 --- a/src/geometries/polytopes.jl +++ b/src/geometries/polytopes.jl @@ -160,10 +160,17 @@ function (c::Chain)(t) if t < 0 || t > 1 throw(DomainError(t, "c(t) is not defined for t outside [0, 1].")) end - n = nvertices(c) - !isclosed(c) - k = min(n - 1, floor(Int, n * t)) - s, _ = iterate(segments(c), k) - s(n * t - k) + segs = segments(c) + measures = cumsum(measure.(segs)) + measures /= measures[end] + # measures[k] <= t < measures[k + 1] + k = searchsortedfirst(measures, t) - 1 + s, _ = iterate(segs, k) + mk = k == 0 ? 0 : measures[k] + mk1 = measures[k + 1] + # Map t to [0, 1] in the segment + t2 = (t - mk) / (mk1 - mk) + s(t2) end # implementations of Chain diff --git a/test/polytopes.jl b/test/polytopes.jl index 9a66e7e1e..7fe876820 100644 --- a/test/polytopes.jl +++ b/test/polytopes.jl @@ -267,16 +267,16 @@ end # parameterization ri = Ring(latlon(45, 0), latlon(45, 90), latlon(90, 90)) @test ri(T(0)) == latlon(45, 0) - @test ri(T(0.25)) == latlon(45, 67.5) - @test ri(T(0.5)) == latlon(67.5, 90) - @test ri(T(0.75)) == latlon(78.75, 67.5) + @test ri(T(0.25)) ≈ latlon(45, 56.25) atol=1e7 * eps(T) * u"m" + @test ri(T(0.5)) ≈ latlon(60, 90) atol=1e7 * eps(T) * u"m" + @test ri(T(0.75)) ≈ latlon(82.5, 75) atol=1e7 * eps(T) * u"m" @test ri(T(1)) == latlon(45, 0) ro = Rope(latlon(45, 0), latlon(45, 90), latlon(90, 90)) @test ro(T(0)) == latlon(45, 0) - @test ro(T(0.25)) == latlon(45, 45) - @test ro(T(0.5)) == latlon(45, 90) - @test ro(T(0.75)) == latlon(67.5, 90) + @test ro(T(0.25)) ≈ latlon(45, 39.375) atol=1e7 * eps(T) * u"m" + @test ro(T(0.5)) ≈ latlon(45, 78.75) atol=1e7 * eps(T) * u"m" + @test ro(T(0.75)) ≈ latlon(63.75, 90) atol=1e7 * eps(T) * u"m" @test ro(T(1)) == latlon(90, 90) ri = Ring(cart.([(1, 1), (2, 2), (3, 3)])) From a45df4fc57e569bb04bc81fb13a868735fbd35d8 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Thu, 19 Dec 2024 12:00:03 +0100 Subject: [PATCH 08/13] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/polytopes.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/polytopes.jl b/test/polytopes.jl index 7fe876820..7810657d3 100644 --- a/test/polytopes.jl +++ b/test/polytopes.jl @@ -267,16 +267,16 @@ end # parameterization ri = Ring(latlon(45, 0), latlon(45, 90), latlon(90, 90)) @test ri(T(0)) == latlon(45, 0) - @test ri(T(0.25)) ≈ latlon(45, 56.25) atol=1e7 * eps(T) * u"m" - @test ri(T(0.5)) ≈ latlon(60, 90) atol=1e7 * eps(T) * u"m" - @test ri(T(0.75)) ≈ latlon(82.5, 75) atol=1e7 * eps(T) * u"m" + @test ri(T(0.25)) ≈ latlon(45, 56.25) atol = 1e7 * eps(T) * u"m" + @test ri(T(0.5)) ≈ latlon(60, 90) atol = 1e7 * eps(T) * u"m" + @test ri(T(0.75)) ≈ latlon(82.5, 75) atol = 1e7 * eps(T) * u"m" @test ri(T(1)) == latlon(45, 0) ro = Rope(latlon(45, 0), latlon(45, 90), latlon(90, 90)) @test ro(T(0)) == latlon(45, 0) - @test ro(T(0.25)) ≈ latlon(45, 39.375) atol=1e7 * eps(T) * u"m" - @test ro(T(0.5)) ≈ latlon(45, 78.75) atol=1e7 * eps(T) * u"m" - @test ro(T(0.75)) ≈ latlon(63.75, 90) atol=1e7 * eps(T) * u"m" + @test ro(T(0.25)) ≈ latlon(45, 39.375) atol = 1e7 * eps(T) * u"m" + @test ro(T(0.5)) ≈ latlon(45, 78.75) atol = 1e7 * eps(T) * u"m" + @test ro(T(0.75)) ≈ latlon(63.75, 90) atol = 1e7 * eps(T) * u"m" @test ro(T(1)) == latlon(90, 90) ri = Ring(cart.([(1, 1), (2, 2), (3, 3)])) From ae84e58a8400a6befb8744a14730a324c4cd08b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= <julio.hoffimann@gmail.com> Date: Thu, 19 Dec 2024 09:29:29 -0300 Subject: [PATCH 09/13] Fix type instabilities --- src/geometries/polytopes.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/geometries/polytopes.jl b/src/geometries/polytopes.jl index bd922408e..463011bbc 100644 --- a/src/geometries/polytopes.jl +++ b/src/geometries/polytopes.jl @@ -161,16 +161,16 @@ function (c::Chain)(t) throw(DomainError(t, "c(t) is not defined for t outside [0, 1].")) end segs = segments(c) - measures = cumsum(measure.(segs)) - measures /= measures[end] - # measures[k] <= t < measures[k + 1] - k = searchsortedfirst(measures, t) - 1 + sums = cumsum(measure.(segs)) + sums /= last(sums) + # find k such that sums[k] ≤ t < sums[k + 1] + k = searchsortedfirst(sums, t) - 1 + # select segment s at index k s, _ = iterate(segs, k) - mk = k == 0 ? 0 : measures[k] - mk1 = measures[k + 1] - # Map t to [0, 1] in the segment - t2 = (t - mk) / (mk1 - mk) - s(t2) + # reparametrization of t within s + ∑ₖ = iszero(k) ? zero(eltype(sums)) : sums[k] + ∑ₖ₊₁ = sums[k + 1] + s((t - ∑ₖ) / (∑ₖ₊₁ - ∑ₖ)) end # implementations of Chain From 44c227040cb63b1b0c33df0fc9ff57dc8c434d81 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <joshua.lampert@uni-hamburg.de> Date: Thu, 19 Dec 2024 13:44:00 +0100 Subject: [PATCH 10/13] add tests for Cartesian coords --- test/polytopes.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/test/polytopes.jl b/test/polytopes.jl index 7810657d3..38e63f5c0 100644 --- a/test/polytopes.jl +++ b/test/polytopes.jl @@ -272,6 +272,12 @@ end @test ri(T(0.75)) ≈ latlon(82.5, 75) atol = 1e7 * eps(T) * u"m" @test ri(T(1)) == latlon(45, 0) + ri = boundary(Box(cart(0,0), cart(1,1))) + @test ri(T(0)) == cart(0, 0) + @test ri(T(0.25)) == cart(1, 0) + @test ri(T(0.5)) == cart(1, 1) + @test ri(T(0.75)) == cart(0, 1) + ro = Rope(latlon(45, 0), latlon(45, 90), latlon(90, 90)) @test ro(T(0)) == latlon(45, 0) @test ro(T(0.25)) ≈ latlon(45, 39.375) atol = 1e7 * eps(T) * u"m" @@ -279,6 +285,13 @@ end @test ro(T(0.75)) ≈ latlon(63.75, 90) atol = 1e7 * eps(T) * u"m" @test ro(T(1)) == latlon(90, 90) + ro = Rope(cart(0,0), cart(3, 0), cart(4,0)) + @test ro(T(0)) == cart(0, 0) + @test ro(T(0.25)) == cart(1, 0) + @test ro(T(0.5)) == cart(2, 0) + @test ro(T(0.75)) == cart(3, 0) + @test ro(T(1)) == cart(4, 0) + ri = Ring(cart.([(1, 1), (2, 2), (3, 3)])) ro = Rope(cart.([(1, 1), (2, 2), (3, 3)])) @test sprint(show, ri) == "Ring((x: 1.0 m, y: 1.0 m), (x: 2.0 m, y: 2.0 m), (x: 3.0 m, y: 3.0 m))" From 08c112aa7fa78b8ec38539b8c8bf962d1cb8cda4 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Thu, 19 Dec 2024 13:45:18 +0100 Subject: [PATCH 11/13] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/polytopes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/polytopes.jl b/test/polytopes.jl index 38e63f5c0..cdc403944 100644 --- a/test/polytopes.jl +++ b/test/polytopes.jl @@ -272,7 +272,7 @@ end @test ri(T(0.75)) ≈ latlon(82.5, 75) atol = 1e7 * eps(T) * u"m" @test ri(T(1)) == latlon(45, 0) - ri = boundary(Box(cart(0,0), cart(1,1))) + ri = boundary(Box(cart(0, 0), cart(1, 1))) @test ri(T(0)) == cart(0, 0) @test ri(T(0.25)) == cart(1, 0) @test ri(T(0.5)) == cart(1, 1) @@ -285,7 +285,7 @@ end @test ro(T(0.75)) ≈ latlon(63.75, 90) atol = 1e7 * eps(T) * u"m" @test ro(T(1)) == latlon(90, 90) - ro = Rope(cart(0,0), cart(3, 0), cart(4,0)) + ro = Rope(cart(0, 0), cart(3, 0), cart(4, 0)) @test ro(T(0)) == cart(0, 0) @test ro(T(0.25)) == cart(1, 0) @test ro(T(0.5)) == cart(2, 0) From bcaf0f261475c278e91e8602d77ee761b4386766 Mon Sep 17 00:00:00 2001 From: Joshua Lampert <joshua.lampert@uni-hamburg.de> Date: Thu, 19 Dec 2024 13:46:20 +0100 Subject: [PATCH 12/13] remove approximate latlon tests --- test/polytopes.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/test/polytopes.jl b/test/polytopes.jl index 38e63f5c0..cfe71b205 100644 --- a/test/polytopes.jl +++ b/test/polytopes.jl @@ -267,9 +267,6 @@ end # parameterization ri = Ring(latlon(45, 0), latlon(45, 90), latlon(90, 90)) @test ri(T(0)) == latlon(45, 0) - @test ri(T(0.25)) ≈ latlon(45, 56.25) atol = 1e7 * eps(T) * u"m" - @test ri(T(0.5)) ≈ latlon(60, 90) atol = 1e7 * eps(T) * u"m" - @test ri(T(0.75)) ≈ latlon(82.5, 75) atol = 1e7 * eps(T) * u"m" @test ri(T(1)) == latlon(45, 0) ri = boundary(Box(cart(0,0), cart(1,1))) @@ -280,9 +277,6 @@ end ro = Rope(latlon(45, 0), latlon(45, 90), latlon(90, 90)) @test ro(T(0)) == latlon(45, 0) - @test ro(T(0.25)) ≈ latlon(45, 39.375) atol = 1e7 * eps(T) * u"m" - @test ro(T(0.5)) ≈ latlon(45, 78.75) atol = 1e7 * eps(T) * u"m" - @test ro(T(0.75)) ≈ latlon(63.75, 90) atol = 1e7 * eps(T) * u"m" @test ro(T(1)) == latlon(90, 90) ro = Rope(cart(0,0), cart(3, 0), cart(4,0)) From e4011bdb787d45e8ea8f8908aaa1075d81ba0179 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= <julio.hoffimann@gmail.com> Date: Thu, 19 Dec 2024 09:50:05 -0300 Subject: [PATCH 13/13] Adjust tests --- test/polytopes.jl | 37 +++++++++++++++++-------------------- 1 file changed, 17 insertions(+), 20 deletions(-) diff --git a/test/polytopes.jl b/test/polytopes.jl index 890fe6717..d677e6433 100644 --- a/test/polytopes.jl +++ b/test/polytopes.jl @@ -265,26 +265,23 @@ end @test crs(centroid(r)) === crs(r) # parameterization - ri = Ring(latlon(45, 0), latlon(45, 90), latlon(90, 90)) - @test ri(T(0)) == latlon(45, 0) - @test ri(T(1)) == latlon(45, 0) - - ri = boundary(Box(cart(0, 0), cart(1, 1))) - @test ri(T(0)) == cart(0, 0) - @test ri(T(0.25)) == cart(1, 0) - @test ri(T(0.5)) == cart(1, 1) - @test ri(T(0.75)) == cart(0, 1) - - ro = Rope(latlon(45, 0), latlon(45, 90), latlon(90, 90)) - @test ro(T(0)) == latlon(45, 0) - @test ro(T(1)) == latlon(90, 90) - - ro = Rope(cart(0, 0), cart(3, 0), cart(4, 0)) - @test ro(T(0)) == cart(0, 0) - @test ro(T(0.25)) == cart(1, 0) - @test ro(T(0.5)) == cart(2, 0) - @test ro(T(0.75)) == cart(3, 0) - @test ro(T(1)) == cart(4, 0) + r = boundary(Box(cart(0, 0), cart(1, 1))) + @test r(T(0)) == cart(0, 0) + @test r(T(0.25)) == cart(1, 0) + @test r(T(0.5)) == cart(1, 1) + @test r(T(0.75)) == cart(0, 1) + r = Rope(cart(0, 0), cart(3, 0), cart(4, 0)) + @test r(T(0)) == cart(0, 0) + @test r(T(0.25)) == cart(1, 0) + @test r(T(0.5)) == cart(2, 0) + @test r(T(0.75)) == cart(3, 0) + @test r(T(1)) == cart(4, 0) + r = Ring(latlon(45, 0), latlon(45, 90), latlon(90, 90)) + @test r(T(0)) == latlon(45, 0) + @test r(T(1)) == latlon(45, 0) + r = Rope(latlon(45, 0), latlon(45, 90), latlon(90, 90)) + @test r(T(0)) == latlon(45, 0) + @test r(T(1)) == latlon(90, 90) ri = Ring(cart.([(1, 1), (2, 2), (3, 3)])) ro = Rope(cart.([(1, 1), (2, 2), (3, 3)]))