Small improvements and cleanup and delaunayedges_fast

JuliaGeometry/VoronoiDelaunay.jl#47

Project.toml

@@ -7,5 +7,4 @@ version = "0.1.0"
 GeometricalPredicates = "fd0ad045-b25c-564e-8f9c-8ef5c5f21267"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 VoronoiDelaunay = "72f80fcb-8c52-57d9-aff0-40c1a3526986"

src/GRPF.jl

@@ -10,7 +10,6 @@ to this git repo.
 module GRPF
 
 using LinearAlgebra
-using Statistics
 using StaticArrays
 import GeometricalPredicates: intriangle
 using VoronoiDelaunay
@@ -77,20 +76,18 @@ in the vicinity of a root or pole.
 
 Notes:
  - Order of `𝓔` is not guaranteed.
- - Count of `phasediffs` of value 1 and 3 can differ from Matlab in normal operation,
+ - Count of phasediffs `ΔQ` of value 1 and 3 can differ from Matlab in normal operation,
  because it depends on "direction" of edge.
 """
 function candidateedges(
     tess::DelaunayTessellation2D{IndexablePoint2D},
     quadrants::Vector{Int8}
     )
 
-    # `phasediffs` for diagnosis and plotting only
-    @debug phasediffs = Vector{Int8}()
     𝓔 = Vector{DelaunayEdge{IndexablePoint2D}}()
 
-    for edge in delaunayedges(tess)
-        e::DelaunayEdge{IndexablePoint2D} = edge  # TODO: infer edge type
+    @inbounds for edge in delaunayedges_fast(tess)
+        e = edge::DelaunayEdge{IndexablePoint2D}
         nodea, nodeb = geta(e), getb(e)
         idxa, idxb = getindex(nodea), getindex(nodeb)
 
@@ -100,13 +97,11 @@ function candidateedges(
         #     idxa, idxb = idxb, idxa
         # end
 
-        ΔQ = mod(quadrants[idxa] - quadrants[idxb], 4)
-        ΔQ == 2 && push!(𝓔, e)
-
-        @debug push!(phasediffs, ΔQ)
+        ΔQ = mod(quadrants[idxa] - quadrants[idxb], 4)  # phase difference
+        if ΔQ == 2
+            push!(𝓔, e)
+        end
     end
-
-    @debug return 𝓔, phasediffs
     return 𝓔
 end
 
@@ -172,7 +167,8 @@ function zone1newnodes!(
     n1a = geta(triangle1)
     n1b = getb(triangle1)
     push!(newnodes, (n1a+n1b)/2)
-    for ii = 1:length(triangles)-1
+
+    @inbounds for ii = 1:length(triangles)-1
         na = geta(triangles[ii])
         nb = getb(triangles[ii])
         nc = getc(triangles[ii])
@@ -203,7 +199,9 @@ function addnewnode!(
     if distance(geom2fcn(node1), geom2fcn(node2)) > tolerance
         avgnode = (node1+node2)/2
         for ii in eachindex(newnodes)
-            distance(newnodes[ii], avgnode) < 2eps() && return nothing
+            if distance(newnodes[ii], avgnode) < 2*eps()
+                return nothing
+            end
         end
         push!(newnodes, avgnode)  # only executed if we haven't already returned
     end
@@ -218,8 +216,8 @@ function zone2newnodes!(
     triangles::Vector{DelaunayTriangle{IndexablePoint2D}}
     )
 
-    for ii in eachindex(triangles)
-        triangle = triangles[ii]::DelaunayTriangle{IndexablePoint2D}
+    @inbounds for triangle in triangles
+        # triangle = triangles[ii]::DelaunayTriangle{IndexablePoint2D}
         na = geta(triangle)
         nb = getb(triangle)
         nc = getc(triangle)
@@ -248,18 +246,20 @@ function contouredges(
 
     # Edges of triangles that contain at least 1 of `edges`
     tmpedges = Vector{DelaunayEdge{IndexablePoint2D}}()
-    for triangle in tess
+    @inbounds for triangle in tess
         # We don't know which "direction" the edges are defined in the triangle,
         # so we need to test both
-        edgea = DelaunayEdge(geta(triangle), getb(triangle))
-        edgearev = DelaunayEdge(getb(triangle), geta(triangle))
-        edgeb = DelaunayEdge(getb(triangle), getc(triangle))
-        edgebrev = DelaunayEdge(getc(triangle), getb(triangle))
-        edgec = DelaunayEdge(getc(triangle), geta(triangle))
-        edgecrev = DelaunayEdge(geta(triangle), getc(triangle))
+        pa, pb, pc = geta(triangle), getb(triangle), getc(triangle)
+
+        edgea = DelaunayEdge(pa, pb)
+        edgearev = DelaunayEdge(pb, pa)
+        edgeb = DelaunayEdge(pb, pc)
+        edgebrev = DelaunayEdge(pc, pb)
+        edgec = DelaunayEdge(pc, pa)
+        edgecrev = DelaunayEdge(pa, pc)
 
         # Does triangle contain edge?
-        for edge in edges
+        @inbounds for edge in edges
             if edgea == edge || edgeb == edge || edgec == edge ||
                 edgearev == edge || edgebrev == edge || edgecrev == edge
                 push!(tmpedges, edgea, edgeb, edgec)
@@ -271,7 +271,7 @@ function contouredges(
     # Remove duplicate (reverse) edges from `tmpedges` and otherwise append to `𝐶`
     𝐶 = Vector{DelaunayEdge{IndexablePoint2D}}()
     duplicateedges = zeros(Int, length(tmpedges))
-    for (idxa, edgea) in enumerate(tmpedges)
+    @inbounds for (idxa, edgea) in enumerate(tmpedges)
         if duplicateedges[idxa] == 0
             for (idxb, edgeb) in enumerate(tmpedges)
                 # Check if Edge(a,b) == Edge(b, a), i.e. there are duplicate edges
@@ -304,7 +304,7 @@ function splittriangles(
     zone1triangles = Vector{DelaunayTriangle{IndexablePoint2D}}()
     zone2triangles = Vector{DelaunayTriangle{IndexablePoint2D}}()
     ii = 0
-    for triangle in tess
+    @inbounds for triangle in tess
         ii += 1
         if trianglecounts[ii] > 1
             push!(zone1triangles, triangle)
@@ -369,6 +369,7 @@ function evaluateregions!(
     refnode = getb(𝐶[1])
     popfirst!(𝐶)
     while length(𝐶) > 0
+        # TODO: redesign inside this while loop so nextedgeidx isn't Union{Int, Array{Int}}
         nextedgeidx = findall(e->geta(e)==refnode, 𝐶)
 
         if isempty(nextedgeidx)
@@ -408,30 +409,37 @@ function rootsandpoles(
     geom2fcn::Function
     )
 
-    numregions = size(regions, 1)
-    q = Vector{Union{Missing, Int}}(undef, numregions)
-    z = Vector{ComplexF64}(undef, numregions)
-    for ii in eachindex(regions)
-        # XXX: ORDER OF REGIONS??? XXX
-        quadrantsequence = [convert(Int8, quadrants[getindex(node)]) for node in regions[ii]]
+    numregions = length(regions)
+
+    zroots = Vector{ComplexF64}()
+    zpoles = Vector{ComplexF64}()
+    @inbounds for ii in eachindex(regions)
+        # XXX: Order of regions?
+        quadrantsequence = [quadrants[getindex(node)] for node in regions[ii]]
+
         # Sign flip because `regions[ii]` are in opposite order of Matlab??
         dQ = -diff(quadrantsequence)
-        for jj in eachindex(dQ)
-            dQ[jj] == 3 && (dQ[jj] = -1)
-            dQ[jj] == -3 && (dQ[jj] = 1)
-            # ``|ΔQ| = 2`` is ambiguous; cannot tell whether phase increases or decreases by two quadrants
-            abs(dQ[jj]) == 2 && (dQ[jj] = 0)
+        @inbounds for jj in eachindex(dQ)
+            if dQ[jj] == 3
+                dQ[jj] = -1
+            elseif dQ[jj] == -3
+                dQ[jj] = 1
+            elseif abs(dQ[jj]) == 2
+                # ``|ΔQ| = 2`` is ambiguous; cannot tell whether phase increases or decreases by two quadrants
+                dQ[jj] = 0
+            end
         end
-        q[ii] = sum(dQ)/4
-        z[ii] = mean(geom2fcn.(regions[ii]))
-    end
-    zroots = [z[i] for i in eachindex(z) if q[i] > 0]
-    zroots_multiplicity = filter(x->x>0, q)
+        q = sum(dQ)/4
+        z = sum(geom2fcn.(regions[ii]))/length(regions[ii])
 
-    zpoles = [z[i] for i in eachindex(z) if q[i] < 0]
-    zpoles_multiplicity = filter(x->x<0, q)
+        if q > 0
+            push!(zroots, z)
+        elseif q < 0
+            push!(zpoles, z)
+        end
+    end
 
-    return zroots, zroots_multiplicity, zpoles, zpoles_multiplicity
+    return zroots, zpoles
 end
 
 """
@@ -447,7 +455,7 @@ function tesselate!(
 
     # Initialize
     numnodes = tess._total_points_added
-    # @assert numnodes == 0
+    @assert numnodes == 0
 
     𝓔 = Vector{DelaunayEdge{IndexablePoint2D}}()
     quadrants = Vector{Int8}()
@@ -466,23 +474,15 @@ function tesselate!(
         numnodes += numnewnodes
 
         # Determine candidate edges that may be near a root or pole
-        @debug 𝓔, phasediffs = candidateedges(tess, quadrants)
         𝓔 = candidateedges(tess, quadrants)
         isempty(𝓔) && error("No roots in the domain")
 
         # Select candidate edges that are longer than the chosen tolerance
         select𝓔 = filter(e -> longedge(e, tolerance, geom2fcn), 𝓔)
-        @debug isempty(select𝓔) && return tess, 𝓔, quadrants, phasediffs
         isempty(select𝓔) && return tess, 𝓔, quadrants
 
         max𝓔length = maximum(distance(geom2fcn(e)) for e in select𝓔)
 
-        @debug begin
-            "Candidate edges length max: $max𝓔length"
-            max𝓔length < tolerance && return tess, 𝓔, quadrants, phasediffs
-        end
-        max𝓔length < tolerance && return tess, 𝓔, quadrants
-
         # How many times does each triangle contain a `select𝓔` node?
         trianglecounts = counttriangleswithnodes(tess, select𝓔)
         zone1triangles, zone2triangles = splittriangles(tess, trianglecounts)
@@ -498,7 +498,6 @@ function tesselate!(
         setindex!.(newnodes, (1:length(newnodes)).+numnodes)
     end
 
-    @debug return tess, 𝓔, quadrants, phasediffs
     return tess, 𝓔, quadrants
 end
 
@@ -516,9 +515,9 @@ function grpf(
     tess, 𝓔, quadrants = tesselate!(tess, newnodes, fcn, geom2fcn, tolerance)
     𝐶 = contouredges(tess, 𝓔)
     regions = evaluateregions!(𝐶, geom2fcn)
-    zroots, zroots_multiplicity, zpoles, zpoles_multiplicity = rootsandpoles(regions, quadrants, geom2fcn)
+    zroots, zpoles = rootsandpoles(regions, quadrants, geom2fcn)
 
-    return zroots, zroots_multiplicity, zpoles, zpoles_multiplicity
+    return zroots, zpoles
 end
 
 end # module

src/GeneralFunctions.jl

@@ -16,7 +16,7 @@ end
 Return true if `edge` has length greater than `tolerance`.
 """
 function longedge(edge::DelaunayEdge, tolerance, geom2fcn::Function)
-    return distance(geom2fcn(edge)...) > tolerance  # TODO: splat performance?
+    return distance(geom2fcn(edge)) > tolerance
 end
 
 """
@@ -53,9 +53,6 @@ function rectangulardomain(Zb::Complex, Ze::Complex, Δr)
     # NOTE: Matlab values of `x` differ slightly because of Matlab's float handling and
     # transpose operator.
     # `sum(x)` is much closer between Julia and Matlab if the above line for `x` is:
-    # `matlabx = [x' ((2:2:m) .- 1)'*dx .+ real(Zb)]'`
-
-    # TEMP
     # x = [x' ((2:2:m) .- 1)'*dx .+ real(Zb)]'
     # x = reshape(x, length(x))
 

src/VoronoiDelaunayExtensions.jl

@@ -16,3 +16,28 @@ Base.getindex(p::IndexablePoint2D) = p._index
 Base.setindex!(p::IndexablePoint2D, v::Int) = setfield!(p, :_index, v)
 Base.:+(p1::IndexablePoint2D, p2::IndexablePoint2D) = IndexablePoint2D(getx(p1)+getx(p2), gety(p1)+gety(p2), -1)
 Base.:/(p1::IndexablePoint2D, n::Real) = IndexablePoint2D(getx(p1)/n, gety(p1)/n, -1)
+
+
+# Improved performance of `delaunayedges()`
+# see: https://github.com/JuliaGeometry/VoronoiDelaunay.jl/issues/47
+function delaunayedges_fast(t::DelaunayTessellation2D)
+    result = DelaunayEdge[]
+    for ix in 2:t._last_trig_index
+        tr = t._trigs[ix]
+        isexternal(tr) && continue
+
+        ix_na = tr._neighbour_a
+        if (ix_na > ix) || isexternal(t._trigs[ix_na])
+            push!(result, DelaunayEdge(getb(tr), getc(tr)))
+        end
+        ix_nb = tr._neighbour_b
+        if (ix_nb > ix) || isexternal(t._trigs[ix_nb])
+            push!(result, DelaunayEdge(geta(tr), getc(tr)))
+        end
+        ix_nc = tr._neighbour_c
+        if (ix_nc > ix) || isexternal(t._trigs[ix_nc])
+            push!(result, DelaunayEdge(geta(tr), getb(tr)))
+        end
+    end
+    return result
+end

test/ComplexModes.jl

@@ -44,7 +44,7 @@ function complexmodes(z)
       w = det(W)
 end
 
-R = 1.
+R = 1.0
 r = 0.15
 
 origcoords = diskdomain(R, r)
@@ -69,7 +69,7 @@ tess, 𝓔, quadrants = GRPF.tesselate!(tess, newnodes, pt -> complexmodes(geom2
 𝐶 = GRPF.contouredges(tess, 𝓔)
 regions = GRPF.evaluateregions!(𝐶, e -> geom2fcn(e, ra, rb, ia, ib))
 
-zroots, zroots_multiplicity, zpoles, zpoles_multiplicity = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
+zroots, zpoles = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
 
 sort!(zroots, by = x -> (real(x), imag(x)))
 sort!(zpoles, by = x -> (real(x), imag(x)))

test/DefaultFunction.jl

@@ -11,10 +11,10 @@ function defaultfcn(z)
 end
 
 # Analysis parameters
-xb = -2.  # real part begin
-xe = 2.  # real part end
-yb = -2.  # imag part begin
-ye = 2.  # imag part end
+xb = -2  # real part begin
+xe = 2  # real part end
+yb = -2  # imag part begin
+ye = 2  # imag part end
 r = 0.2  # initial mesh step
 tolerance = 1e-9
 
@@ -39,7 +39,7 @@ tess, 𝓔, quadrants = GRPF.tesselate!(tess, newnodes, pt -> defaultfcn(geom2fc
 𝐶 = GRPF.contouredges(tess, 𝓔)
 regions = GRPF.evaluateregions!(𝐶, e -> geom2fcn(e, ra, rb, ia, ib))
 
-zroots, zroots_multiplicity, zpoles, zpoles_multiplicity = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
+zroots, zpoles = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
 
 sort!(zroots, by = x -> (real(x), imag(x)))
 sort!(zpoles, by = x -> (real(x), imag(x)))

test/GrapheneTransmissionLine.jl

@@ -59,7 +59,7 @@ tess, 𝓔, quadrants = GRPF.tesselate!(tess, newnodes, pt -> graphenefunction(g
 𝐶 = GRPF.contouredges(tess, 𝓔)
 regions = GRPF.evaluateregions!(𝐶, e -> geom2fcn(e, ra, rb, ia, ib))
 
-zroots, zroots_multiplicity, zpoles, zpoles_multiplicity = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
+zroots, zpoles = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
 
 sort!(zroots, by = x -> (real(x), imag(x)))
 sort!(zpoles, by = x -> (real(x), imag(x)))

test/LossyMultilayeredWaveguide.jl

@@ -18,10 +18,10 @@ function wvgd(z)
 end
 
 # Analysis parameters
-xb = 1.  # real part begin
+xb = 1.0  # real part begin
 xe = 2.5  # real part end
-yb = -1.  # imag part begin
-ye = 1.  # imag part end
+yb = -1.0  # imag part begin
+ye = 1.0  # imag part end
 r = 0.5  # initial mesh step
 tolerance = 1e-9
 
@@ -46,7 +46,7 @@ tess, 𝓔, quadrants = GRPF.tesselate!(tess, newnodes, pt -> wvgd(geom2fcn(pt,
 𝐶 = GRPF.contouredges(tess, 𝓔)
 regions = GRPF.evaluateregions!(𝐶, e -> geom2fcn(e, ra, rb, ia, ib))
 
-zroots, zroots_multiplicity, zpoles, zpoles_multiplicity = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
+zroots, zpoles = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
 
 sort!(zroots, by = x -> (real(x), imag(x)))
 sort!(zpoles, by = x -> (real(x), imag(x)))

test/SimpleRationalFunction.jl

@@ -3,10 +3,10 @@ function simplefcn(z)
 end
 
 # Analysis parameters
-xb = -2.  # real part begin
-xe = 2.  # real part end
-yb = -2.  # imag part begin
-ye = 2.  # imag part end
+xb = -2  # real part begin
+xe = 2  # real part end
+yb = -2  # imag part begin
+ye = 2  # imag part end
 r = 0.1  # initial mesh step
 tolerance = 1e-9
 
@@ -31,7 +31,7 @@ tess, 𝓔, quadrants = GRPF.tesselate!(tess, newnodes, pt -> simplefcn(geom2fcn
 𝐶 = GRPF.contouredges(tess, 𝓔)
 regions = GRPF.evaluateregions!(𝐶, e -> geom2fcn(e, ra, rb, ia, ib))
 
-zroots, zroots_multiplicity, zpoles, zpoles_multiplicity = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
+zroots, zpoles = GRPF.rootsandpoles(regions, quadrants, e -> geom2fcn(e, ra, rb, ia, ib))
 
 sort!(zroots, by = x -> (real(x), imag(x)))