augustt198 / Quickhull.jl

MIT License
15 stars 2 forks source link

Fix type instability, speeding up the code a lot #4

Closed IlianPihlajamaa closed 3 weeks ago

IlianPihlajamaa commented 3 weeks ago

For 10^4 points, timing the facets call with and without the change:

function new_delaunay_facets(dhull)
    hull = dhull.hull
    maxlift = maximum(last, hull.pts)

    return filter(hull.facets.arr) do facet
        D = length(first(hull.pts))
        above_pt = sum(hull.pts[i] for i in facet.plane.point_indices) / D
        above_pt = setindex(above_pt, above_pt[end] + 2maxlift, D)
        Quickhull.hyperplane_dist(facet.plane, above_pt, hull.pts) < 0
    end
end

function new_facets(dhull)
    return Quickhull.mappedarray(f -> Quickhull.NgonFace(f.plane.point_indices), new_delaunay_facets(dhull))
end

using BenchmarkTools
r = rand(2, 10000)
r2 = rand(SVector{2, Float64}, 10000)

tri = Quickhull.delaunay(r)
tri2 = Quickhull.delaunay(r2)

println("Benchmarking facets new")
facs_new = @btime new_facets($tri)
facs2_new = @btime new_facets($tri2)

println("Benchmarking facets old")
facs_old = @btime Quickhull.facets($tri)
facs2_old = @btime Quickhull.facets($tri2)

@assert facs_new == facs_old
@assert facs2_new == facs2_old

Benchmarking facets new 769.400 μs (2 allocations: 156.30 KiB) 390.100 μs (2 allocations: 156.30 KiB) Benchmarking facets old 76.063 ms (977732 allocations: 39.18 MiB) 74.385 ms (977782 allocations: 38.87 MiB)

augustt198 commented 3 weeks ago

Good catch, thanks! I'll have to look into the 2x performance difference still