Minimize allocation of single element arrays

[garborg]

Only matters for the functions that do less work, but cuts allocation quite a bit, and good to apply consistently.

[yeesian]

oo, I haven't seen this style of allocating arrays before. do you have any benchmarks? I'm okay with merging it either way!

[garborg]

julia> module Test

       using LibGEOS
       using LibGEOS: GEOSGeom, GEOSDistance, GEOSGeomGetY

       function getGeomY1(ptr::GEOSGeom)
           y = Array(Float64, 1)
           result = GEOSGeomGetY(ptr, pointer(y))
           if result == -1
               error("LibGEOS: Error in GEOSGeomGetY")
           end
           y[1]
       end

       let y = Array(Float64, 1)
           global getGeomY2
           function getGeomY2(ptr::GEOSGeom)
               result = GEOSGeomGetY(ptr, pointer(y))
               if result == -1
                   error("LibGEOS: Error in GEOSGeomGetY")
               end
               y[1]
           end
       end

       let y = Array(Float64, 1), ptrd = pointer(y)
           global getGeomY3
julia> module Test

       using LibGEOS
       using LibGEOS: GEOSGeom, GEOSDistance, GEOSGeomGetY

       function getGeomY1(ptr::GEOSGeom)
           y = Array(Float64, 1)
           result = GEOSGeomGetY(ptr, pointer(y))
           if result == -1
               error("LibGEOS: Error in GEOSGeomGetY")
           end
           y[1]
       end

       let y = Array(Float64, 1)
           global getGeomY2
           function getGeomY2(ptr::GEOSGeom)
               result = GEOSGeomGetY(ptr, pointer(y))
               if result == -1
                   error("LibGEOS: Error in GEOSGeomGetY")
               end
               y[1]
           end
       end

       let y = Array(Float64, 1), ptrd = pointer(y)
           global getGeomY3
           function getGeomY3(ptr::GEOSGeom)
               result = GEOSGeomGetY(ptr, ptrd)
               if result == -1
                   error("LibGEOS: Error in GEOSGeomGetY")
               end
               y[1]
           end
       end

       function test_d1(n)
           x = 0.0
           for _ = 1:n
               x += getGeomY1(pt1)
           end
           x
       end

       function test_d2(n)
           x = 0.0
           for _ = 1:n
               x += getGeomY2(pt1)
           end
           x
       end

       function test_d3(n)
           x = 0.0
           for _ = 1:n
               x += getGeomY3(pt1)
           end
           x
       end

       pt1 = LibGEOS.createPoint(1.1, 2.2)
       pt2 = LibGEOS.createPoint(2.1, 3.2)

       test_d1(2)
       test_d2(2)
       test_d3(2)

       n = 10_000_000
       for _ = 1:3
           println()
           gc(); @time test_d1(n)
           gc(); @time test_d2(n)
           gc(); @time test_d3(n)
       end

       end
Warning: replacing module Test

   1.176 seconds      (20000 k allocations: 916 MB, 9.78% gc time)
 666.456 milliseconds (10000 k allocations: 153 MB, 1.62% gc time)
 672.659 milliseconds (10000 k allocations: 153 MB, 1.68% gc time)

   1.148 seconds      (20000 k allocations: 916 MB, 9.33% gc time)
 640.996 milliseconds (10000 k allocations: 153 MB, 1.72% gc time)
 654.454 milliseconds (10000 k allocations: 153 MB, 1.69% gc time)

   1.136 seconds      (20000 k allocations: 916 MB, 9.45% gc time)
 673.100 milliseconds (10000 k allocations: 153 MB, 1.73% gc time)
 681.891 milliseconds (10000 k allocations: 153 MB, 1.70% gc time)

For something that does more work, like distance, I only saw a ~6% speed up even though the allocation stats are similar.

[garborg]

OT: But, I haven't looked at how the package works with abstract points. Is it free to pass in a geometry made up of LL{T <: Datum), or (Float64, Float64, Float64)?

I'm curious because ideally, we'd want to avoid using Array{Float64,1} for points except in cases where we just haven't customized deserialization from a standard format yet. I don't know how much work it will be to get there, though.

[yeesian]

OT: But, I haven't looked at how the package works with abstract points. Is it free to pass in a geometry made up of LL{T <: Datum), or (Float64, Float64, Float64)?

That was meant to be the role of https://github.com/JuliaGeo/GeoInterface.jl/blob/master/src/GeoInterface.jl#L21-L37. The idea was for LLA/ENU/etc to implement AbstractPosition, and have it work with everything else. I haven't yet begun on that work. And we'll still need to construct LibGEOS types from GeoInterface types though.

I'm curious because ideally, we'd want to avoid using Array{Float64,1} for points except in cases where we just haven't customized deserialization from a standard format yet. I don't know how much work it will be to get there, though.

Agreed. I'm talking to @sjkelly and @SimonDanisch about JuliaGeometry/GeometryTypes.jl#3.

[yeesian]

I'm happy with the changes in this PR. Shall we merge it?

[garborg]

Great, let's merge it.

I'm glad you're out there hashing out how how this corner of the ecosystem is going to look! Wish I could have made it.

I have an old branch of Geodesy.jl that abused the type system a bit to get points to carry their coordinate system for free (minimizing method lookups and no memory overhead). I'll push it up and send you a comment -- not sure it'd fit into the work you're doing with @sjkelly and @SimonDanisch, but it was certainly nice for calculating distances, converting between coordinate systems, etc., without adding extra method parameters up and down through the call stack.

[yeesian]

I have an old branch of Geodesy.jl that abused the type system a bit to get points to carry their coordinate system for free (minimizing method lookups and no memory overhead).

I'll be really interested to see it. I had also spent some thinking cycles about map projections, and (if i find time) provide wrappers for OGR (see https://groups.google.com/d/msg/julia-geo/oPz2_gz9K-g/LVfou9f93AEJ).

[garborg]

That would be great to get into Julia.

I just pushed the branch up in JuliaGeo/Geodesy.jl#11 -- not wedded to anything specific about the design, except that having points parameterized on their coordinate system somehow seems very useful.