Faster hash and equality (#159)

* overload hash and equality

* Update src/geos_types.jl

* bump version

* wip

* fix hash

* implement == via equals

* implement coordinatewise equality

* overload Base.isapprox

* more isapprox tests

* somewhat faster compare

* faster

* speedup hash and equality

* test that hash eq don't allocate much

* eliminate more allocs

* add Random as test dependency

* fix

* fix

* relax alloc tests

* use MersenneTwister instead of Xoshiro in tests

* simplify

* fix ci

* add another early return to compare

* fix

* fix

* add some coordinates! tests

* fix typo

* get rid of Random test dependency

Author: jw3126

src/geos_functions.jl

@@ -292,22 +292,31 @@ function getZ(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())
     zcoords
 end
 
-function getCoordinates(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = get_global_context())
+function getCoordinates!(out::Vector{Float64}, ptr::GEOSCoordSeq, i::Integer, context)
     if !(0 < i <= getSize(ptr, context))
         error(
             "LibGEOS: i=$i is out of bounds for CoordSeq with size=$(getSize(ptr, context))",
         )
     end
     ndim = getDimensions(ptr, context)
-    coord = Array{Float64}(undef, ndim)
-    start = pointer(coord)
-    floatsize = sizeof(Float64)
-    GEOSCoordSeq_getX_r(context, ptr, i - 1, start)
-    GEOSCoordSeq_getY_r(context, ptr, i - 1, start + floatsize)
-    if ndim == 3
-        GEOSCoordSeq_getZ_r(context, ptr, i - 1, start + 2 * floatsize)
+    @assert length(out) == ndim
+    GC.@preserve out begin
+        start = pointer(out)
+        floatsize = sizeof(Float64)
+        GEOSCoordSeq_getX_r(context, ptr, i - 1, start)
+        GEOSCoordSeq_getY_r(context, ptr, i - 1, start + floatsize)
+        if ndim == 3
+            GEOSCoordSeq_getZ_r(context, ptr, i - 1, start + 2 * floatsize)
+        end
     end
-    coord
+    out
+end
+
+function getCoordinates(ptr::GEOSCoordSeq, i::Integer, context::GEOSContext = get_global_context())
+    ndim = getDimensions(ptr, context)
+    out = Array{Float64}(undef, ndim)
+    getCoordinates!(out, ptr, i, context)
+    out
 end
 
 function getCoordinates(ptr::GEOSCoordSeq, context::GEOSContext = get_global_context())

src/geos_types.jl

@@ -89,7 +89,7 @@ mutable struct LineString <: AbstractGeometry
         finalizer(destroyGeom, line)
         line
     end
-    #create a linestring from a list of coordiantes
+    #create a linestring from a list of coordinates
     function LineString(coords::Vector{Vector{Float64}}, context::GEOSContext = get_global_context())
         line = new(createLineString(coords, context), context)
         finalizer(destroyGeom, line)
@@ -116,7 +116,7 @@ mutable struct MultiLineString <: AbstractGeometry
         finalizer(destroyGeom, multiline)
         multiline
     end
-    # create a multilinestring from a list of linestring coordiantes
+    # create a multilinestring from a list of linestring coordinates
     MultiLineString(multiline::Vector{Vector{Vector{Float64}}},context::GEOSContext = get_global_context()) =
         MultiLineString(
             createCollection(
@@ -309,6 +309,31 @@ const geomtypes = [
     GeometryCollection,
 ]
 
+const HasCoordSeq = Union{LineString, LinearRing, Point}
+"""
+
+    coordinates!(out::Vector{Float64}, geo::$HasCoordSeq, i::Integer)
+    coordinates!(out::Vector{Float64}, geo::Point)
+
+Copy the coordinates of the ith point of geo into `out`.
+"""
+function coordinates!(out, geo::HasCoordSeq, i::Integer, ctx::GEOSContext=get_context(geo))
+    GC.@preserve out geo begin
+        seq = GEOSGeom_getCoordSeq_r(ctx, geo)::Ptr
+        getCoordinates!(out, seq, i, ctx)
+    end
+end
+function coordinates!(out, geo::Point, ctx::GEOSContext=get_context(geo))
+    coordinates!(out, geo, 1, ctx)
+end
+
+function has_coord_seq(geo::AbstractGeometry)
+    false
+end
+function has_coord_seq(::HasCoordSeq)
+    true
+end
+
 Base.@kwdef struct IsApprox
     atol::Float64=0.0
     rtol::Float64=sqrt(eps(Float64))
@@ -323,51 +348,78 @@ end
 function Base.isapprox(geo1::AbstractGeometry, geo2::AbstractGeometry; kw...)::Bool
     compare(IsApprox(;kw...), geo1, geo2)
 end
-function (cmp::IsApprox)(geo1::AbstractGeometry, geo2::AbstractGeometry)::Bool
-    compare(cmp, geo1, geo2)
-end
-function compare(≅, geo1::AbstractGeometry, geo2::AbstractGeometry)::Bool
+function compare(cmp, geo1::AbstractGeometry, geo2::AbstractGeometry, ctx=get_context(geo1))::Bool
     (typeof(geo1) === typeof(geo2)) || return false
-    ng1 = ngeom(geo1)
-    ng2 = ngeom(geo2)
-    ng1 == ng2 || return false
-    for i in 1:ng1
-        (getgeom(geo1, i) ≅ getgeom(geo2, i)) || return false
+    if (geo1 === geo2) && (cmp === isequal)
+        return true
     end
-    return true
-end
-function compare(≅, pt1::Point, pt2::Point)::Bool
-    is3d = GeoInterface.is3d(pt1)
-    is3d === GeoInterface.is3d(pt2) || return false
-    GeoInterface.getcoord(pt1,1) ≅ GeoInterface.getcoord(pt2,1) || return false
-    GeoInterface.getcoord(pt1,2) ≅ GeoInterface.getcoord(pt2,2) || return false
-    if is3d
-        GeoInterface.getcoord(pt1,3) ≅ GeoInterface.getcoord(pt2,3) || return false
+    if has_coord_seq(geo1)
+        return compare_coord_seqs(cmp, geo1, geo2, ctx)
+    else
+        ng1 = ngeom(geo1)
+        ng2 = ngeom(geo2)
+        ng1 == ng2 || return false
+        for i in 1:ng1
+            compare(cmp, getgeom(geo1, i), getgeom(geo2, i), ctx) || return false
+        end
     end
     return true
 end
 
-function _norm(x,y,z)
-    return sqrt(x^2 + y^2 + z^2)
+npoints(pt::Point) = 1
+npoints(geo::HasCoordSeq) = GeoInterface.ngeom(geo)
+
+function compare_coord_seqs(cmp, geo1, geo2, ctx)
+    ncoords1 = GeoInterface.ncoord(geo1)
+    ncoords2 = GeoInterface.ncoord(geo2)
+    ncoords1 == ncoords2 || return false
+    ncoords1 == 0 && return true
+    np1 = npoints(geo1)
+    np2 = npoints(geo2)
+    np1 == np2 || return false
+    coords1 = Vector{Float64}(undef, ncoords1)
+    coords2 = Vector{Float64}(undef, ncoords1)
+    for i in 1:np1
+        coordinates!(coords1, geo1, i, ctx)
+        coordinates!(coords2, geo2, i, ctx)
+        cmp(coords1, coords2) || return false
+    end
+    return true
 end
 
-function compare(cmp::IsApprox, pt1::Point, pt2::Point)::Bool
-    is3d = GeoInterface.is3d(pt1)
-    is3d === GeoInterface.is3d(pt2) || return false
-    x1 = GeoInterface.getcoord(pt1,1)
-    x2 = GeoInterface.getcoord(pt2,1)
-    y1 = GeoInterface.getcoord(pt1,2)
-    y2 = GeoInterface.getcoord(pt2,2)
-    if is3d
-        z1 = GeoInterface.getcoord(pt1,3) 
-        z2 = GeoInterface.getcoord(pt2,3)
-    else
-        z1 = 0.0
-        z2 = 0.0
+function compare_coord_seqs(cmp::IsApprox, geo1, geo2, ctx)
+    ncoords1 = GeoInterface.ncoord(geo1)
+    ncoords2 = GeoInterface.ncoord(geo2)
+    ncoords1 == ncoords2 || return false
+    ncoords1 == 0 && return true
+    @assert ncoords1 in 2:3
+    ncoords1 == ncoords2 || return false
+    np1 = npoints(geo1)
+    np2 = npoints(geo2)
+    np1 == np2 || return false
+    coords1 = Vector{Float64}(undef, ncoords1)
+    coords2 = Vector{Float64}(undef, ncoords1)
+    # isapprox of two vectors is calculated using their euclidean norms and the norm of their difference
+    # we compute (the squares) of these norms in an allocation free way
+    s1 = 0.0
+    s2 = 0.0
+    s12 = 0.0
+    for i in 1:np1
+        coordinates!(coords1, geo1, i, ctx)
+        coordinates!(coords2, geo2, i, ctx)
+        if ncoords1 == 2
+            x1,y1 = coords1
+            x2,y2 = coords2
+            s12 += (x1 - x2)^2 + (y1 - y2)^2
+        else
+            x1,y1,z1 = coords1
+            x2,y2,z2 = coords2
+            s12 += (x1 - x2)^2 + (y1 - y2)^2 + (z1 - z2)^2
+        end
+        s1 += sum(abs2, coords1)
+        s2 += sum(abs2, coords2)
     end
-    lhs = _norm(x1 - x2, y1 - y2, z1 - z2)
-    rhs = cmp.atol + cmp.rtol * max(_norm(x1,y1,z2), _norm(x2,y2,z2))
-    return lhs <= rhs
+    return sqrt(s12) <= cmp.atol + cmp.rtol * sqrt(max(s1, s2))
 end
 
 typesalt(::Type{GeometryCollection} ) = 0xd1fd7c6403c36e5b
@@ -381,25 +433,31 @@ typesalt(::Type{Point}              ) = 0x4b5c101d3843160e
 typesalt(::Type{Polygon}            ) = 0xa5c895d62ef56723
 
 function Base.hash(geo::AbstractGeometry, h::UInt)::UInt
-    salt = typesalt(typeof(geo))
-    h = hash(salt, h)
-    for i in 1:ngeom(geo)
-        g = getgeom(geo,i)
-        h = hash(g, h)
+    h = hash(typesalt(typeof(geo)), h)
+    if has_coord_seq(geo) 
+        return hash_coord_seq(geo, h)
+    else
+        for i in 1:ngeom(geo)
+            h = hash(getgeom(geo, i), h)
+        end
     end
     return h
 end
-
-function Base.hash(pt::Point, h::UInt)::UInt
-    h = hash(getcoord(pt,1), h)
-    h = hash(getcoord(pt,2), h)
-    if GeoInterface.is3d(pt)
-        h = hash(getcoord(pt,3), h)
+function hash_coord_seq(geo::HasCoordSeq, h::UInt)::UInt
+    nc = GeoInterface.ncoord(geo)
+    if nc == 0
+        return h
+    end
+    buf = Vector{Float64}(undef, nc)
+    ctx = get_context(geo)
+    for i in 1:npoints(geo)
+        coordinates!(buf, geo, i, ctx)
+        h = hash(buf, h)
     end
-    h = hash(getcoord(pt,2), h)
     return h
 end
 
+
 # teach ccall how to get the pointer to pass to libgeos
 # this way the Julia compiler will track the lifetimes for us
 Base.unsafe_convert(::Type{Ptr{Cvoid}}, x::AbstractGeometry) = x.ptr

test/test_misc.jl

@@ -46,6 +46,17 @@ end
     @test LibGEOS.intersects(p2, q1, ctx3)
 end
 
+@testset "coordinates!" begin
+    coordinates! = LibGEOS.coordinates!
+    buf3 = zeros(3)
+    buf2 = zeros(2)
+    @test coordinates!(buf3, readgeom("POINT(1 2 3)")) == Float64[1,2,3] == buf3
+    @test coordinates!(buf2, readgeom("POINT(1 2)")) == Float64[1,2] == buf2
+    @test coordinates!(buf2, readgeom("POINT(1 2)"), 1) == Float64[1,2] == buf2
+    @test coordinates!(buf2, readgeom("LINESTRING (130 240, 650 240)"), 1) == buf2 == [130, 240]
+    @test coordinates!(buf2, readgeom("LINESTRING (130 240, 650 240)"), 2) == buf2 == [650, 240]
+end
+
 @testset "hash eq" begin
     ctx1 = LibGEOS.GEOSContext()
     ctx2 = LibGEOS.GEOSContext()
@@ -66,6 +77,15 @@ end
     @test readgeom("LINESTRING (130 240, 650 240)") ≈ readgeom("LINESTRING (130 240, -650 240)") atol=1300
     @test readgeom("LINESTRING (130 240, 650 240)") ≈ readgeom("LINESTRING (130 240, 650 240.00000001)")
 
+    @test isapprox(readgeom("POLYGON((1 1,0 0,            1 2,2 2,2 4,1 1))"),
+                   readgeom("POLYGON((1 1,0 0.00000000001,1 2,2 2,2 4,1 1))"))
+
+    pt = readgeom("POINT(-1 NaN)")
+    @test isequal(pt, pt)
+    @test pt != pt
+    @test !(isapprox(pt, pt))
+    @test !(isapprox(pt, pt, atol=Inf))
+
     pt = readgeom("POINT(0 NaN)")
     @test isequal(pt, pt)
     @test pt != pt
@@ -115,6 +135,40 @@ end
     end
 end
 
+@testset "performance hash eq" begin
+    pts1 = [[sin(x), cos(x), 1] for x in range(0, 2pi, length=1000)]
+    pts1[end] = pts1[1]
+    lr1 = LinearRing(pts1)
+    pts2 = copy(pts1)
+    pts2[453] = 2*pts2[453]
+    lr2 = LinearRing(pts2)
+    @test !(lr1 == lr2)
+    @test !isequal(lr1, lr2)
+    @test !isapprox(lr1, lr2)
+    
+    @test 300 > @allocated lr1 == lr2
+    @test 300 > @allocated isequal(lr1, lr2)
+    @test 300 > @allocated isapprox(lr1, lr2)
+
+    hash(lr1) != hash(lr2)
+    @test 300 > @allocated hash(lr1)
+
+    poly1 = Polygon(lr1)
+    poly2 = Polygon(lr2)
+
+    poly2 = LinearRing(pts2)
+    @test !(poly1 == poly2)
+    @test !isequal(poly1, poly2)
+    @test !isapprox(poly1, poly2)
+    
+    @test 300 > @allocated poly1 == poly2
+    @test 300 > @allocated isequal(poly1, poly2)
+    @test 300 > @allocated isapprox(poly1, poly2)
+
+    @test hash(poly1) != hash(poly2)
+    @test 300 > @allocated hash(poly1)
+end
+
 @testset "show it like you build it" begin
     for geo in [
         readgeom("POINT(0 0)")