Fix copy and deepcopy for PolygonTrees and PolygonHierarchys

NEWS.md

@@ -1,14 +1,15 @@
 # Changelog
 
 ## 1.6.0
-- Feature: Define `reverse` for `AbstractParametricCurve`s, making it easier to reverse the orientation of a curve. See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
-- Bugfix: Fixed an issue with `LineSegment` not returning the exact endpoints at `t=1`, which can be problematic when joining boundary nodes. This has been fixed. See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
-- Bugfix: Introduced `is_linear` to fix issues with boundary enrichment of domains with `LineSegment`s. In particular, `LineSegment`s are no longer enriched. See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
-- Bugfix: `orientation_markers` now uses `uniquetol` instead of `unique` for the final set of markers (it already did it for the intermediate markers). See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
-- Bugfix: For large `Tuple`s, functions like `eval_fnc_at_het_tuple_two_elements` are problematic and allocate more than their non-type-stable counterparts. To get around this, for `Tuple`s of length `N > 32`, the non-type-stable version is used. See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
-- Bigfix: Fixed issue with `use_barriers` when a ghost edge is selected at random during point location. See [#196](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/196).
-- Feature: Introduced the (currently internal) function `get_positive_curve_indices` for finding curves with positive orientation in a `Triangulation`. [#196](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/196).
+- Define `reverse` for `AbstractParametricCurve`s, making it easier to reverse the orientation of a curve. See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
+- Fixed an issue with `LineSegment` not returning the exact endpoints at `t=1`, which can be problematic when joining boundary nodes. This has been fixed. See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
+- Introduced `is_linear` to fix issues with boundary enrichment of domains with `LineSegment`s. In particular, `LineSegment`s are no longer enriched. See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
+- `orientation_markers` now uses `uniquetol` instead of `unique` for the final set of markers (it already did it for the intermediate markers). See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
+- For large `Tuple`s, functions like `eval_fnc_at_het_tuple_two_elements` are problematic and allocate more than their non-type-stable counterparts. To get around this, for `Tuple`s of length `N > 32`, the non-type-stable version is used. See [#195](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/195).
+- Fixed issue with `use_barriers` when a ghost edge is selected at random during point location. See [#196](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/196).
+- Introduced the (currently internal) function `get_positive_curve_indices` for finding curves with positive orientation in a `Triangulation`. [#196](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/196).
 - `is_exterior_curve`, `is_interior_curve`, `num_exterior_curves`, and `is_disjoint` are now defined based on `get_positive_curve_indices` rather than `get_exterior_curve_indices`. See [#196](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/196).
+- `copy` and `deepcopy` are now correctly implemented for `PolygonTree`s and `PolygonHierarchy`s. See [#199](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/199)
 
 ## 1.5.0
 

src/data_structures/trees/polygon_hierarchy.jl

@@ -16,19 +16,19 @@
 Constructs a [`PolygonTree`](@ref) with `parent`, `index`, and `height`, and no children.
 """
 mutable struct PolygonTree{I}
-    parent::Union{Nothing, PolygonTree{I}}
+    parent::Union{Nothing,PolygonTree{I}}
     @const children::Set{PolygonTree{I}}
     @const index::I
     height::Int
 end
-PolygonTree{I}(parent::Union{Nothing, PolygonTree{I}}, index, height) where {I} = PolygonTree{I}(parent, Set{PolygonTree{I}}(), index, height)
+PolygonTree{I}(parent::Union{Nothing,PolygonTree{I}}, index, height) where {I} = PolygonTree{I}(parent, Set{PolygonTree{I}}(), index, height)
 function hash_tree(tree::PolygonTree)
     height = get_height(tree)
     index = get_index(tree)
     parent_index = has_parent(tree) ? get_index(get_parent(tree)) : 0
     h = hash((parent_index, index, height))
     children = collect(get_children(tree))
-    sort!(children, by = get_index)
+    sort!(children, by=get_index)
     for child in children
         h = hash((h, hash_tree(child)))
     end
@@ -46,22 +46,23 @@
     height1 ≠ height2 && return false
     return hash_tree(tree1) == hash_tree(tree2)
 end
-@static if VERSION ≥ v"1.10"
-    function Base.deepcopy(tree::PolygonTree) # without this definition, deepcopy would occassionally segfault 
-        @warn "Deepcopy on PolygonTrees is not currently supported. Returning the tree without copying." maxlog = 1 # https://github.com/JuliaGeometry/DelaunayTriangulation.jl/issues/129
-        return tree
-        #=
-        parent = get_parent(tree)
-        children = get_children(tree)
-        index = get_index(tree)
-        height = get_height(tree)
-        new_parent = isnothing(parent) ? nothing : deepcopy(parent)
-        new_children = deepcopy(children)
-        new_index = deepcopy(index)
-        return PolygonTree(new_parent, new_children, new_index, height)
-        =# # Why does the above still segfault sometimes?
-        # TODO: Fix
+function Base.deepcopy_internal(tree::PolygonTree{I}, dict::IdDict) where {I}
+    haskey(dict, tree) && return dict[tree]
+    copy_tree = PolygonTree{I}(nothing, get_index(tree), get_height(tree))
+    dict[tree] = copy_tree
+    new_parent = Base.deepcopy_internal(get_parent(tree), dict)
+    if !isnothing(new_parent)
+        set_parent!(copy_tree, new_parent)
     end
+    for child in get_children(tree)
+        add_child!(copy_tree, Base.deepcopy_internal(child, dict))
+    end
+    return copy_tree
+end
+function Base.copy(tree::PolygonTree{I}) where {I}
+    parent = get_parent(tree)
+    new_parent = isnothing(parent) ? nothing : copy(parent)
+    return PolygonTree{I}(new_parent, copy(get_children(tree)), copy(get_index(tree)), copy(get_height(tree)))
 end
 function Base.show(io::IO, ::MIME"text/plain", tree::PolygonTree)
     print(io, "PolygonTree at height $(get_height(tree)) with index $(get_index(tree)) and $(length(get_children(tree))) children")
@@ -158,7 +159,7 @@
 - `polygon_orientations::BitVector`: A `BitVector` of length `n` where `n` is the number of polygons in the hierarchy. The `i`th entry is `true` if the `i`th polygon is positively oriented, and `false` otherwise.
 - `bounding_boxes::Vector{BoundingBox}`: A `Vector` of [`BoundingBox`](@ref)s of length `n` where `n` is the number of polygons in the hierarchy. The `i`th entry is the [`BoundingBox`](@ref) of the `i`th polygon.
 - `trees::Dict{I,PolygonTree{I}}`: A `Dict` mapping the index of a polygon to its [`PolygonTree`](@ref). The keys of `trees` are the roots of each individual tree, i.e. the outer-most polygons.
-- `reorder_cache::Vector{PolygonTree{I}}`: A `Vector used for caching trees to be deleted in [`reorder_subtree!`](@ref).
+- `reorder_cache::Vector{PolygonTree{I}}`: A `Vector` used for caching trees to be deleted in [`reorder_subtree!`](@ref).
 
 !!! note "One-based indexing"
 
@@ -174,29 +175,40 @@
 struct PolygonHierarchy{I}
     polygon_orientations::BitVector
     bounding_boxes::Vector{BoundingBox}
-    trees::Dict{I, PolygonTree{I}}
+    trees::Dict{I,PolygonTree{I}}
     reorder_cache::Vector{PolygonTree{I}}
 end
-PolygonHierarchy{I}() where {I} = PolygonHierarchy{I}(BitVector(), BoundingBox[], Dict{I, PolygonTree{I}}(), PolygonTree{I}[])
-@static if VERSION ≥ v"1.10"
-    function Base.deepcopy(hierarchy::PolygonHierarchy{I}) where {I} # without this definition, deepcopy would occassionally segfault 
-        polygon_orientations = get_polygon_orientations(hierarchy)
-        bounding_boxes = get_bounding_boxes(hierarchy)
-        trees = get_trees(hierarchy)
-        reorder_cache = get_reorder_cache(hierarchy)
-        new_polygon_orientations = copy(polygon_orientations)
-        new_bounding_boxes = copy(bounding_boxes)
-        new_trees = Dict{I, PolygonTree{I}}()
-        for (index, tree) in trees
-            new_trees[index] = deepcopy(tree)
-        end
-        new_reorder_cache = similar(reorder_cache)
-        for (index, tree) in enumerate(reorder_cache)
-            new_reorder_cache[index] = deepcopy(tree)
-        end
-        return PolygonHierarchy{I}(new_polygon_orientations, new_bounding_boxes, new_trees, new_reorder_cache)
+PolygonHierarchy{I}() where {I} = PolygonHierarchy{I}(BitVector(), BoundingBox[], Dict{I,PolygonTree{I}}(), PolygonTree{I}[])
+function Base.deepcopy_internal(hierarchy::PolygonHierarchy{I}, dict::IdDict) where {I}
+    haskey(dict, hierarchy) && return dict[hierarchy]
+    copy_polygon_orientations = Base.deepcopy_internal(get_polygon_orientations(hierarchy), dict)
+    copy_bounding_boxes = Base.deepcopy_internal(get_bounding_boxes(hierarchy), dict)
+    copy_trees = Dict{I, PolygonTree{I}}()
+    for (key, tree) in get_trees(hierarchy)
+        copy_trees[key] = Base.deepcopy_internal(tree, dict)
+    end
+    copy_reorder_cache = Vector{PolygonTree{I}}()
+    for tree in get_reorder_cache(hierarchy)
+        push!(copy_reorder_cache, Base.deepcopy_internal(tree, dict))
     end
+    copy_hierarchy = PolygonHierarchy(
+        copy_polygon_orientations,
+        copy_bounding_boxes,
+        copy_trees,
+        copy_reorder_cache
+    )
+    dict[hierarchy] = copy_hierarchy
+    return copy_hierarchy
 end
+function Base.copy(hierarchy::PolygonHierarchy{I}) where {I}
+    return PolygonHierarchy(
+        copy(get_polygon_orientations(hierarchy)),
+        copy(get_bounding_boxes(hierarchy)),
+        copy(get_trees(hierarchy)),
+        copy(get_reorder_cache(hierarchy))
+    )
+end
+
 function Base.empty!(hierarchy::PolygonHierarchy)
     polygon_orientations = get_polygon_orientations(hierarchy)
     bounding_boxes = get_bounding_boxes(hierarchy)
@@ -288,7 +300,7 @@
 Returns the indices of the positively oriented curves of `hierarchy` as a generator, i.e. 
 as a lazy result.
 """
-function get_positive_curve_indices(hierarchy::PolygonHierarchy) 
+function get_positive_curve_indices(hierarchy::PolygonHierarchy)
     orientations = get_polygon_orientations(hierarchy)
     return (index for (index, orientation) in enumerate(orientations) if orientation)
 end
@@ -365,7 +377,7 @@
 
 Returns a [`PolygonHierarchy`](@ref) defining the polygon hierarchy for a given set of `points`. This defines a hierarchy with a single polygon.
 """
-function construct_polygon_hierarchy(points; IntegerType = Int)
+function construct_polygon_hierarchy(points; IntegerType=Int)
     hierarchy = PolygonHierarchy{IntegerType}()
     return construct_polygon_hierarchy!(hierarchy, points)
 end
@@ -385,11 +397,11 @@
 Returns a [`PolygonHierarchy`](@ref) defining the polygon hierarchy for a given set of `boundary_nodes` that define a set of piecewise 
 linear curves. 
 """
-function construct_polygon_hierarchy(points, boundary_nodes; IntegerType = Int)
+function construct_polygon_hierarchy(points, boundary_nodes; IntegerType=Int)
     hierarchy = PolygonHierarchy{IntegerType}()
     return construct_polygon_hierarchy!(hierarchy, points, boundary_nodes)
 end
-construct_polygon_hierarchy(points, ::Nothing; IntegerType = Int) = construct_polygon_hierarchy(points; IntegerType)
+construct_polygon_hierarchy(points, ::Nothing; IntegerType=Int) = construct_polygon_hierarchy(points; IntegerType)
 function construct_polygon_hierarchy!(hierarchy::PolygonHierarchy{I}, points, boundary_nodes) where {I}
     if !has_boundary_nodes(boundary_nodes)
         return construct_polygon_hierarchy!(hierarchy, points)
@@ -611,7 +623,7 @@
 
 Expands the bounding boxes of `hierarchy` by a factor of `perc` in each direction.
 """
-function expand_bounds!(hierarchy::PolygonHierarchy, perc = 0.1)
+function expand_bounds!(hierarchy::PolygonHierarchy, perc=0.1)
     bboxes = get_bounding_boxes(hierarchy)
     for (i, bbox) in enumerate(bboxes)
         bboxes[i] = expand(bbox, perc)
@@ -634,9 +646,9 @@
 - `IntegerType=Int`: The integer type to use for indexing the polygons.
 - `n=4096`: The number of points to use for filling in the boundary curves in [`polygonise`](@ref).
 """
-function construct_polygon_hierarchy(points, boundary_nodes, boundary_curves; IntegerType = Int, n = 4096)
+function construct_polygon_hierarchy(points, boundary_nodes, boundary_curves; IntegerType=Int, n=4096)
     new_points, new_boundary_nodes = polygonise(points, boundary_nodes, boundary_curves; n)
     hierarchy = PolygonHierarchy{IntegerType}()
     return construct_polygon_hierarchy!(hierarchy, new_points, new_boundary_nodes)
 end
-construct_polygon_hierarchy(points, boundary_nodes, ::Tuple{}; IntegerType = Int, n = 4096) = construct_polygon_hierarchy(points, boundary_nodes; IntegerType)
+construct_polygon_hierarchy(points, boundary_nodes, ::Tuple{}; IntegerType=Int, n=4096) = construct_polygon_hierarchy(points, boundary_nodes; IntegerType)

test/data_structures/polygon_hierarchy.jl

@@ -42,7 +42,7 @@ for _ in 1:20 # Run many times to make sure the segfault is gone
     ]
 
     curve_IV = [CircularArc((1.0, 0.0), (1.0, 0.0), (0.0, 0.0))]
-    points_IV = NTuple{2, Float64}[]
+    points_IV = NTuple{2,Float64}[]
 
     curve_V = [BezierCurve([(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)])]
     points_V = [(0.0, 0.0), (0.2, 0.25)]
@@ -73,13 +73,13 @@ for _ in 1:20 # Run many times to make sure the segfault is gone
 
     curve_IX =
         [
-        [
-            [1, 2, 3, 4, 5, 6, 7, 1],
-        ],
-        [
-            [CircularArc((0.6, 0.5), (0.6, 0.5), (0.5, 0.5), positive = false)],
-        ],
-    ]
+            [
+                [1, 2, 3, 4, 5, 6, 7, 1],
+            ],
+            [
+                [CircularArc((0.6, 0.5), (0.6, 0.5), (0.5, 0.5), positive=false)],
+            ],
+        ]
     points_IX = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.5, 1.5), (0.0, 1.0), (0.0, 0.5), (0.0, 0.2)]
 
     curve_X = [
@@ -111,7 +111,7 @@ for _ in 1:20 # Run many times to make sure the segfault is gone
             [12, 11, 10, 12],
         ],
         [
-            [CircularArc((1.1, -3.0), (1.1, -3.0), (0.0, -3.0), positive = false)],
+            [CircularArc((1.1, -3.0), (1.1, -3.0), (0.0, -3.0), positive=false)],
         ],
     ]
     points_XI = [(-2.0, 0.0), (0.0, 0.0), (2.0, 0.0), (-2.0, -5.0), (2.0, -5.0), (2.0, -1 / 10), (-2.0, -1 / 10), (-1.0, -3.0), (0.0, -4.0), (0.0, -2.3), (-0.5, -3.5), (0.9, -3.0)]
@@ -335,11 +335,16 @@ for _ in 1:20 # Run many times to make sure the segfault is gone
         @test traverse_tree(hierarchy.trees[7]) ≠ traverse_tree(hierarchy.trees[15])
         @test compare_trees(hierarchy, hierarchy)
         @test compare_trees(deepcopy(hierarchy), deepcopy(hierarchy))
+        @test compare_trees(hierarchy, deepcopy(hierarchy))
+        @test compare_trees(hierarchy, copy(hierarchy))
+        @test compare_trees(copy(hierarchy), copy(hierarchy))
+        @test !(hierarchy === deepcopy(hierarchy))
         @test !compare_trees(hierarchy, hierarchy2)
         @test hierarchy.trees[7] ≠ hierarchy.trees[15]
         @test hierarchy == hierarchy
         @test deepcopy(hierarchy) == hierarchy
         @test deepcopy(hierarchy) == deepcopy(hierarchy)
+        @test !(hierarchy === deepcopy(hierarchy))
         @test hierarchy ≠ hierarchy2
         @test deepcopy(hierarchy) ≠ hierarchy2
     end
@@ -357,7 +362,7 @@ for _ in 1:20 # Run many times to make sure the segfault is gone
         hierarchy = DT.construct_polygon_hierarchy(points, nnew_boundary_nodes, boundary_curves)
         DT.expand_bounds!(hierarchy, DT.ε(Float64))
         @test DT.get_bounding_boxes(hierarchy) ⊢ [DT.BoundingBox(-1 - 2DT.ε(Float64), 1 + 2DT.ε(Float64), -1 - 2DT.ε(Float64), 1 + 2DT.ε(Float64))] &&
-            DT.get_polygon_orientations(hierarchy) ⊢ BitVector([1])
+              DT.get_polygon_orientations(hierarchy) ⊢ BitVector([1])
         trees = DT.get_trees(hierarchy)
         @test length(trees) == 1
         tree = trees[1]
@@ -459,4 +464,4 @@ for _ in 1:20 # Run many times to make sure the segfault is gone
         tree2 = trees[2]
         @test DT.get_index(tree2) == 2 && DT.get_height(tree2) == 0 && !DT.has_parent(tree2) && !DT.has_children(tree2)
     end
-end
+end