Extents are somehow very slow

[asinghvi17]

On Julia 1.12

Setup

using Extents, GADM, Chairmarks
import GeometryOps as GO

switzerland = GI.geometry(GI.getfeature(GADM.get("CHE"), 1)) |> GO.tuples

Result

julia> @be GO.prepare_naturally($(switzerland.geom[1]))
Benchmark: 17 samples with 1 evaluation
 min    74.916 μs (29 allocs: 1.349 MiB)
 median 143.375 μs (29 allocs: 1.349 MiB)
 mean   12.350 ms (29.24 allocs: 1.349 MiB, 11.73% gc time)
 max    149.128 ms (33 allocs: 1.349 MiB, 99.86% gc time)

julia> @be GI.extent($(switzerland.geom[1].geom[1].geom |> GI.LinearRing); fallback = true)
Benchmark: 416 samples with 1 evaluation
 min    200.292 μs
 median 207.834 μs
 mean   225.360 μs
 max    384.791 μs

This means that creating a natural tree which allocates etc. is 1.3x faster than gettting the extent which is pretty insane. We should try to fix this somehow, because that means that any geometry without cached extent is getting clobbered in benchmarks.

[rafaqz]

Sounds like you need to profile;)

Very likely the new julia version has introduced type instability?

[asinghvi17]

With type instability I would expect allocations though, but yeah will have to profile and see what's going on.

I think the mapreduce changes made it into 1.12, that may be a part of it

[evetion]

Seems to be mapreduce yes.

edit: But I see the same performance on 1.11 though

[asinghvi17]

Ah it calls getpoint directly, it's not recursive. Maybe that's the speed difference? It also goes over the iterator thrice, once per dimension, which is...not the greatest of ideas (IMO)

GeoInterface.jl/src/fallbacks.jl

Line 134 in 5a85ddc

points = getpoint(t, geom)

[asinghvi17]

Ok so this is a bit long BUT it gets runtime down by nearly 8x while allocating a bit more

julia> @be GI.extent($(switzerland); fallback = true)
Benchmark: 310 samples with 1 evaluation
 min    287.791 μs (15 allocs: 496 bytes)
 median 304.834 μs (15 allocs: 496 bytes)
 mean   306.004 μs (15 allocs: 496 bytes)
 max    395.917 μs (15 allocs: 496 bytes)

julia> @be find_extent($(switzerland))
Benchmark: 2139 samples with 1 evaluation
 min    36.583 μs (25 allocs: 992 bytes)
 median 39.334 μs (25 allocs: 992 bytes)
 mean   39.563 μs (25 allocs: 992 bytes)
 max    255.167 μs (25 allocs: 992 bytes)

import GeometryOps as GO, GeoInterface as GI

function _find_extent(trait::GI.AbstractGeometryTrait, geom, hasZ::GO.BoolsAsTypes)
    points = GI.getpoint(trait, geom)
    xmin, ymin = 0.0, 0.0
    xmax, ymax = 0.0, 0.0
    if GO.istrue(hasZ)
        zmin, zmax = 0.0, 0.0
    end
    for point in points
        x, y = GI.x(point), GI.y(point)
        if x < xmin
            xmin = x
        end
        if x > xmax
            xmax = x
        end
        if y < ymin
            ymin = y
        end
        if y > ymax
            ymax = y
        end
        if GO.istrue(hasZ)
            z = GI.z(point)
            if z < zmin
                zmin = z
            end
            if z > zmax
                zmax = z
            end
        end
    end
    
    if GO.istrue(hasZ)
        return Extent(X = (xmin, xmax), Y = (ymin, ymax), Z = (zmin, zmax))
    else
        return Extent(X=(xmin, xmax), Y=(ymin, ymax))
    end
end
function _find_extent(::GI.PointTrait, geom, hasZ::GO.BoolsAsTypes)
    x = GI.x(geom)
    y = GI.y(geom)
    if GO.istrue(hasZ)
        z = GI.z(geom)
    end
    return GO.istrue(hasZ) ? Extent(X=(x, x), Y=(y, y), Z=(z, z)) : Extent(X=(x, x), Y=(y, y))
end


function find_extent(trait::GI.AbstractTrait, geom, hasZ::GO.BoolsAsTypes)
    f(t, g) = _find_extent(t, g, hasZ)
    return GO.applyreduce(
        GO.WithTrait(f), 
        GO.Extents.union, 
        GO.TraitTarget{Union{GI.AbstractCurveTrait, GI.PointTrait}}(), 
        geom; 
        threaded = false
    )
end

find_extent(trait::GI.AbstractTrait, geom) = find_extent(trait, geom, GO.booltype(GO._is3d(trait, geom)))
find_extent(geom) = find_extent(GI.trait(geom), geom)