Prune implicit machinery flagged by the substrate audit

src/methods/geom_relations/relateng/edge_intersector.jl

@@ -180,15 +180,15 @@ recorded on `computer`.
 - `NestedLoop`: a plain double loop over
   string pairs and segment pairs, with a per-pair segment-extent
   disjointness skip (on `Planar`).
-- Any tree-backed accelerator (e.g. `DoubleSTRtree`): a spatial index
-  (`_relate_edge_index`, currently a `NaturalIndex`) is built over the
-  per-segment extents of each side and traversed with
+- Any tree-backed accelerator (canonically `DoubleNaturalTree`): a spatial
+  index (`_relate_edge_index`, currently a `NaturalIndex`) is built over
+  the per-segment extents of each side and traversed with
   `SpatialTreeInterface.dual_depth_first_search`
   under the `Extents.intersects` predicate.
 - [`AutoAccelerator`](@ref): picks `NestedLoop` below the clipping size
   threshold (`GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS`) and on
-  non-`Planar` manifolds (planar extent trees are not valid there), and the
-  tree path otherwise.
+  manifolds without a segment-extent kernel (neither `Planar` nor
+  `Spherical`), and the tree path otherwise.
 
 After each processed pair `is_result_known(computer)` is consulted and the
 enumeration stops early once the predicate value is determined (the port of
@@ -220,29 +220,20 @@ function process_edge_intersections!(tc::TopologyComputer,
         _select_edge_set_accelerator(m, ssa_list, ssb_list); m, exact)
 end
 
-# STRtrees over planar extents are only valid on the Planar manifold, and
-# below the clipping threshold the nested loop wins anyway.
-function _select_edge_set_accelerator(::Planar, ssa_list, ssb_list)
+# Below the clipping threshold the nested loop wins; above it, the tree
+# path. Valid on both manifolds with a segment-extent kernel: planar boxes
+# on `Planar`, 3D great-circle arc extents (`_segment_extent`) on
+# `Spherical` — `NaturalIndex`, the dual DFS, and `Extents.intersects` are
+# dimension-generic. Other manifolds have no extent kernel and always take
+# the nested loop.
+function _select_edge_set_accelerator(::Union{Planar, Spherical}, ssa_list, ssb_list)
     na = _total_segment_count(ssa_list)
     nb = _total_segment_count(ssb_list)
     if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS &&
             nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
         return NestedLoop()
     else
-        return DoubleSTRtree()
-    end
-end
-# The same threshold heuristic on the sphere: the tree path is valid because
-# the segment extents are 3D great-circle arc extents (`_segment_extent`), and
-# `NaturalIndex` / the dual DFS / `Extents.intersects` are dimension-generic.
-function _select_edge_set_accelerator(::Spherical, ssa_list, ssb_list)
-    na = _total_segment_count(ssa_list)
-    nb = _total_segment_count(ssb_list)
-    if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS &&
-            nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
-        return NestedLoop()
-    else
-        return DoubleSTRtree()
+        return DoubleNaturalTree()
     end
 end
 _select_edge_set_accelerator(::Manifold, ssa_list, ssb_list) = NestedLoop()
@@ -300,7 +291,7 @@ _segment_envs_disjoint(::Manifold, a0, a1, b0, b1) = false
 _relate_edge_index(extents::Vector{<:Extents.Extent}) =
     NaturalIndex(extents; nodecapacity = 16)
 
-# Tree path (any other accelerator, canonically DoubleSTRtree): a spatial
+# Tree path (any other accelerator, canonically DoubleNaturalTree): a spatial
 # index over the per-segment extents of each side, traversed simultaneously.
 function process_edge_intersections!(tc::TopologyComputer,
         ssa_list::AbstractVector{<:RelateSegmentString},

src/methods/geom_relations/relateng/indexed_point_in_area.jl

@@ -48,20 +48,15 @@ Port of JTS `SortedPackedIntervalRTree`, with two representation changes
   of level `k` (an unpaired trailing node is carried up unchanged, as in
   `buildLevel`). The last level is the root.
 - JTS always sorts the leaves by interval midpoint (`NodeComparator`)
-  before packing. Here `sort_leaves = false` skips that and keeps insertion
-  (ring) order — the `NaturalIndexing` observation. Query results are
-  order-independent (every visit is extent-checked), but the layouts trade
-  off differently: midpoint order groups same-`y` segments so a point query
-  descends few subtrees, while a long coastline ring in natural order
-  recrosses the query `y` in many separated runs. Measured on Natural Earth
-  10m Canada: the sort is ~4× the rest of the build, and natural-order
-  queries are ~3× slower. So prepared mode sorts (build once, query
-  forever) and the lazily indexed unprepared path doesn't (its query count
-  is at most a few hundred, far below the ~1000-query crossover; see
-  `locate_on_polygonal`).
+  before packing; so does this port. The sort earns its cost in the
+  index's only (prepared, build-once-query-forever) use: midpoint order
+  groups same-`y` segments so a point query descends few subtrees, where
+  insertion (ring) order recrosses the query `y` in many separated runs —
+  measured on Natural Earth 10m Canada, the sort is ~4× the rest of the
+  build and ring-order queries are ~3× slower.
 """
 struct SortedPackedIntervalRTree{I}
-    # leaf items: midpoint-sorted (`sort_leaves = true`) or insertion order
+    # leaf items, midpoint-sorted
     items::Vector{I}
     # level_min[1][i] / level_max[1][i] is the interval of leaf item i;
     # level k > 1 holds the pairwise-combined extents of level k - 1
@@ -70,19 +65,15 @@ struct SortedPackedIntervalRTree{I}
 end
 
 # Port of insert + init/buildRoot/buildTree/buildLevel, packed eagerly.
-# Without `sort_leaves` the leaf arrays are taken over by the tree, not
-# copied.
 function SortedPackedIntervalRTree(mins::Vector{Float64}, maxs::Vector{Float64},
-        items::Vector{I}; sort_leaves::Bool = true) where {I}
-    if sort_leaves
-        #-- sort the leaf nodes (IntervalRTreeNode.NodeComparator: by
-        #-- midpoint; sortperm is stable, matching Collections.sort)
-        n = length(items)
-        perm = sortperm(Float64[(mins[i] + maxs[i]) / 2 for i in 1:n])
-        mins = mins[perm]
-        maxs = maxs[perm]
-        items = items[perm]
-    end
+        items::Vector{I}) where {I}
+    #-- sort the leaf nodes (IntervalRTreeNode.NodeComparator: by
+    #-- midpoint; sortperm is stable, matching Collections.sort)
+    n = length(items)
+    perm = sortperm(Float64[(mins[i] + maxs[i]) / 2 for i in 1:n])
+    mins = mins[perm]
+    maxs = maxs[perm]
+    items = items[perm]
     level_min = [mins]
     level_max = [maxs]
     #-- now group nodes into blocks of two and build tree up recursively
@@ -285,14 +276,14 @@ struct IndexedPointInAreaLocator{M <: Manifold, E}
     is_empty::Bool
 end
 
-function IndexedPointInAreaLocator(m::Manifold, geom; exact, sort_leaves::Bool = true)
+function IndexedPointInAreaLocator(m::Manifold, geom; exact)
     mins = Float64[]
     maxs = Float64[]
     segs = _PIASegment[]
     n = GI.npoint(geom)
     sizehint!(mins, n); sizehint!(maxs, n); sizehint!(segs, n)
     _interval_index_add_geom!(mins, maxs, segs, GI.trait(geom), geom)
-    index = SortedPackedIntervalRTree(mins, maxs, segs; sort_leaves)
+    index = SortedPackedIntervalRTree(mins, maxs, segs)
     #-- IntervalIndexedGeometry.isEmpty: a (recursively) empty polygonal
     #-- geometry contributes no rings, hence no segments
     return IndexedPointInAreaLocator(m, exact, index, isempty(segs))

src/methods/geom_relations/relateng/point_locator.jl

@@ -235,12 +235,9 @@ bnRule)`; the manifold/`exact` parameters are the only additions (consistent
 with [`AdjacentEdgeLocator`](@ref)). As in JTS, prepared mode swaps the
 per-polygon `SimplePointInAreaLocator` ring loop for a cached
 [`IndexedPointInAreaLocator`](@ref) (indexed_point_in_area.jl), created
-lazily on the first use per polygonal element (Task 22). Unprepared mode
-deviates from Java (which keys indexing on `isPrepared` alone): the first
-query on a polygonal element uses the direct ring loop, but repeat queries
-build and reuse the indexed locator — one O(n) scan beats an O(n) index
-build, while the many area-vertex locations of a multi-element relate
-amortize the index (see `locate_on_polygonal`).
+lazily on the first use per polygonal element (Task 22); unprepared mode
+scans the rings directly on every query. Repeated point location against
+one geometry is what [`prepare`](@ref) is for.
 """
 mutable struct RelatePointLocator{M <: Manifold, E, G, BR <: BoundaryNodeRule, P}
     const m::M
@@ -257,14 +254,10 @@ mutable struct RelatePointLocator{M <: Manifold, E, G, BR <: BoundaryNodeRule, P
     const polygons::Vector{Any}
     const line_boundary::LinearBoundary{BR, P}
     const is_empty::Bool
-    # per-polygonal-element indexed locators, created lazily by
-    # `_get_poly_locator` (Java: polyLocator, filled by getLocator).
-    # Prepared mode fills an entry on its first query; unprepared mode on
-    # its second (see `locate_on_polygonal`).
+    # per-polygonal-element indexed locators (prepared mode only), created
+    # lazily by `_get_poly_locator` on the element's first query
+    # (Java: polyLocator, filled by getLocator)
     const poly_locator::Vector{Union{Nothing, IndexedPointInAreaLocator{M, E}}}
-    # unprepared mode: queries seen per polygonal element, driving the lazy
-    # index heuristic above
-    const poly_query_count::Vector{Int32}
     # lazily built on the first multi-boundary point (Java: adjEdgeLocator)
     adj_edge_locator::Union{Nothing, AdjacentEdgeLocator{M, E, P}}
 end
@@ -285,19 +278,18 @@ function RelatePointLocator(m::Manifold, geom; exact,
     # LinearBoundary behaves identically (no boundary, no boundary points),
     # so it is built unconditionally here.
     line_boundary = LinearBoundary(m, lines, boundary_rule)
-    # Java allocates `polyLocator` for both modes (Simple/Indexed); here both
-    # modes may cache indexed locator objects (unprepared lazily, on repeat
-    # queries), so it is allocated unconditionally.
+    # Java allocates `polyLocator` for both modes (its unprepared arm caches
+    # SimplePointInAreaLocator objects); the direct ring scan here is
+    # stateless, so only prepared mode fills it.
     poly_locator = Vector{Union{Nothing, IndexedPointInAreaLocator{typeof(m), typeof(exact)}}}(
         nothing, length(polygons))
-    poly_query_count = zeros(Int32, length(polygons))
     #-- P cannot be inferred from the `nothing` adj_edge_locator, so spell out
     #-- every type parameter
     return RelatePointLocator{typeof(m), typeof(exact), typeof(geom),
             typeof(boundary_rule), P}(
         m, exact, geom, is_prepared, boundary_rule,
         points, lines, polygons, line_boundary, is_empty, poly_locator,
-        poly_query_count, nothing)
+        nothing)
 end
 
 has_boundary(loc::RelatePointLocator) = has_boundary(loc.line_boundary)
@@ -501,24 +493,9 @@ function locate_on_polygons(loc::RelatePointLocator, p, is_node::Bool, parent_po
     return LOC_EXTERIOR
 end
 
-# Queries a polygonal element absorbs via the direct ring loop before its
-# IndexedPointInAreaLocator is built. Both costs scale with the element's
-# segment count, so one threshold fits all sizes: an unsorted index build
-# costs ~10-13 ring scans (measured on Natural Earth coastlines), making
-# the worst-case regret of switching at 8 about one build. Real relates are
-# bimodal — a handful of queries (barely-touching neighbors, where indexing
-# never pays) or hundreds (one area-vertex location per polygon element of
-# the other geometry), so the threshold rarely sits near the break-even.
-const _LAZY_INDEX_QUERY_THRESHOLD = Int32(8)
-
 # Port of RelatePointLocator.locateOnPolygonal: Java dispatches to a
 # per-polygonal PointOnGeometryLocator — a cached IndexedPointInAreaLocator
-# when prepared, a SimplePointInAreaLocator otherwise. Prepared mode does
-# the same here (Task 22). Unprepared mode deviates from Java: the first
-# query on an element uses the direct SimplePointInAreaLocator ring loop
-# (one O(n) scan beats an O(n) index build + query), but repeat queries —
-# e.g. one area-vertex location per polygon element of the other geometry
-# in a multipolygon/multipolygon relate — build and amortize the index.
+# when prepared, a SimplePointInAreaLocator otherwise. Same here (Task 22).
 function locate_on_polygonal(loc::RelatePointLocator, p, is_node::Bool, parent_polygonal, index::Int)
     polygonal = loc.polygons[index]
     if is_node && parent_polygonal === polygonal
@@ -527,29 +504,18 @@ function locate_on_polygonal(loc::RelatePointLocator, p, is_node::Bool, parent_p
     #-- the RayCrossingCounter horizontal-ray sweep is coordinate-plane
     #-- logic (as is all of JTS), so a future non-planar kernel falls
     #-- through to its own rk_point_in_ring even when prepared
-    if loc.m isa Planar
-        use_index = loc.is_prepared
-        if !use_index
-            count = (loc.poly_query_count[index] += Int32(1))
-            use_index = count > _LAZY_INDEX_QUERY_THRESHOLD
-        end
-        if use_index
-            return locate(_get_poly_locator(loc, index), p)
-        end
+    if loc.is_prepared && loc.m isa Planar
+        return locate(_get_poly_locator(loc, index), p)
     end
     return _locate_point_in_polygonal(loc.m, p, GI.trait(polygonal), polygonal; exact = loc.exact)
 end
 
 # Port of RelatePointLocator.getLocator (indexed arm): lazily create and
-# cache the indexed locator for polygonal element `index`. Prepared mode
-# pays for the midpoint-sorted layout (build once, query forever); the
-# unprepared lazy index skips the sort, which dominates the build cost
-# (see `SortedPackedIntervalRTree`).
+# cache the indexed locator for polygonal element `index`.
 function _get_poly_locator(loc::RelatePointLocator, index::Int)
     locator = loc.poly_locator[index]
     if locator === nothing
-        locator = IndexedPointInAreaLocator(loc.m, loc.polygons[index];
-            exact = loc.exact, sort_leaves = loc.is_prepared)
+        locator = IndexedPointInAreaLocator(loc.m, loc.polygons[index]; exact = loc.exact)
         loc.poly_locator[index] = locator
     end
     return locator

src/methods/geom_relations/relateng/relate_ng.jl

@@ -709,21 +709,15 @@ relate_predicate(p::PreparedRelate, predicate::TopologyPredicate, b) =
 # Whether to prebuild the A-side segment tree, mirroring the dispatch of
 # `process_edge_intersections!` + `_select_edge_set_accelerator`: an explicit
 # `NestedLoop` accelerator never uses a tree; `AutoAccelerator` uses one on
-# `Planar` above the clipping size threshold only (B is unknown at prepare
-# time, so the decision is made on A's segment count alone); any other
-# explicit accelerator always takes the tree path.
+# the manifolds with a segment-extent kernel (`Planar`/`Spherical`) above
+# the clipping size threshold only (B is unknown at prepare time, so the
+# decision is made on A's segment count alone); any other explicit
+# accelerator always takes the tree path.
 _build_prepared_edge_index(m::Manifold, ::IntersectionAccelerator, segs_a) =
     _make_prepared_edge_index(m, segs_a)
 _build_prepared_edge_index(::Manifold, ::NestedLoop, segs_a) = nothing
 _build_prepared_edge_index(::Manifold, ::AutoAccelerator, segs_a) = nothing
-function _build_prepared_edge_index(m::Planar, ::AutoAccelerator, segs_a)
-    _total_segment_count(segs_a) >= GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS ||
-        return nothing
-    return _make_prepared_edge_index(m, segs_a)
-end
-#-- the Spherical tree path is valid too (3D arc extents); above the threshold
-#-- prepared mode indexes A just like Planar
-function _build_prepared_edge_index(m::Spherical, ::AutoAccelerator, segs_a)
+function _build_prepared_edge_index(m::Union{Planar, Spherical}, ::AutoAccelerator, segs_a)
     _total_segment_count(segs_a) >= GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS ||
         return nothing
     return _make_prepared_edge_index(m, segs_a)

src/methods/geom_relations/relateng/topology_predicate.jl

@@ -211,10 +211,9 @@ intersects_exterior_of(p::IMPredicate, is_a::Bool) = is_a ?
     (is_intersects_entry(p, LOC_INTERIOR, LOC_EXTERIOR) || is_intersects_entry(p, LOC_BOUNDARY, LOC_EXTERIOR))
 
 is_intersects_entry(p::IMPredicate, locA, locB) = p.im[locA, locB] >= DIM_P
-# NOTE: unused; kept for JTS IMPredicate API parity. As ported it can never
-# return `false`: matrix entries are initialized to DIM_FALSE and only ever
-# increase, so they never hold DIM_UNKNOWN (-3). Do not use as a real check.
-is_known_entry(p::IMPredicate, locA, locB) = p.im[locA, locB] != DIM_UNKNOWN
+# JTS's isKnownEntry is not ported: entries here are initialized to
+# DIM_FALSE and only ever increase, so they never hold DIM_UNKNOWN and the
+# check could never return false.
 is_dimension_entry(p::IMPredicate, locA, locB, dim) = p.im[locA, locB] == dim
 get_dimension(p::IMPredicate, locA, locB) = p.im[locA, locB]
 

test/methods/relateng/edge_intersector.jl

@@ -346,7 +346,7 @@ ngon(cx, cy, r, n) = GI.Polygon([[
     ]
     for (ga, gb) in fixtures
         tc_truth, pred_truth = run_all_pairs(ga, gb)
-        for acc in (GO.NestedLoop(), GO.DoubleSTRtree(), GO.AutoAccelerator())
+        for acc in (GO.NestedLoop(), GO.DoubleNaturalTree(), GO.AutoAccelerator())
             tc, pred = run_enum(ga, gb, acc)
             @test node_counts(tc) == node_counts(tc_truth)
             @test imstr(pred) == imstr(pred_truth)
@@ -395,7 +395,7 @@ end
     # stop long before enumerating every extent-overlapping pair
     ga = ngon(0.0, 0.0, 1.0, 64)
     gb = ngon(0.1, 0.0, 1.0, 64)
-    for acc in (GO.NestedLoop(), GO.DoubleSTRtree())
+    for acc in (GO.NestedLoop(), GO.DoubleNaturalTree())
         # baseline: the full-matrix predicate is never known early, so its
         # count is the total number of pairs the enumeration would process
         full = run_counted(ga, gb, acc, GO.RelateMatrixPredicate())
@@ -436,7 +436,7 @@ end
     end
 
     counts_loop, im_loop = run_self(GO.NestedLoop())
-    counts_tree, im_tree = run_self(GO.DoubleSTRtree())
+    counts_tree, im_tree = run_self(GO.DoubleNaturalTree())
     @test counts_loop == counts_tree
     @test im_loop == im_tree
     # sanity: the return segment properly crosses all 40 zigzag segments
@@ -461,7 +461,7 @@ end
         GO.process_self_intersections!(tc, ss_list, acc)
         return pred
     end
-    for acc in (GO.NestedLoop(), GO.DoubleSTRtree())
+    for acc in (GO.NestedLoop(), GO.DoubleNaturalTree())
         # baseline: the full-matrix predicate is never known early, so its
         # count is the total number of pairs the enumeration would process
         full = run_self_counted(GO.RelateMatrixPredicate(), acc)