Thread-safety tools with ThreadFunctors, simplified apply via Applicators (#288)

Changelog:
- Refactor the `apply` pipeline to get rid of closures in favour of `Applicator`s, which are able to run `apply` on arrays, geoms and featurecollections.  Some examples are `ApplyToPoint`, `ApplyToGeom`, etc.  These are mainly used by internal methods in the `apply` pipeline.
- Add a `TaskFunctors` struct that allows one to use individual functors per task.  This allows thread safety for reentrant C APIs.
- Refactor `reproject`, in the Proj extension, to use this functionality


* add ThreadFunctors for threaded stateful functors in apply

* include threading.jl

* import ThreadFunctors

* add Applicators

* comments

* Move applicators out of apply.jl

* Replace the apply on pointtrait to polygon guts

with a custom applicator.  Maybe we need a formalized interface for applicators.

* include applicators in core

* add a few more docs in threading.jl

* define applicator interface + remove overwrite

* fix reproject

* add back crs info to rebuilt geoms

* remove error

* add GFT

* add and use applicators

* more applicators

* stop segfault by assigning context post cloning / construction

* fix maptasks on thread functors

* Fix tests

ProjString needs +type=crs but regular String does not.

This is likely because I stopped converting everything to String and started relying on the Proj constructors.

* fixc

* union proj crs

* if statement for stability

* passing tests

* to -> with

---------

Co-authored-by: Anshul Singhvi <[email protected]>

Author: rafaqz

GeometryOpsCore/src/GeometryOpsCore.jl

@@ -27,6 +27,8 @@ include("types/operation.jl")
 include("types/exceptions.jl")
 include("types/booltypes.jl")
 include("types/traittarget.jl")
+include("types/applicators.jl")
+include("types/tasks.jl")
 
 include("apply.jl")
 include("applyreduce.jl")

GeometryOpsCore/src/apply.jl

@@ -145,19 +145,16 @@ end
     _apply(f, target, GI.trait(geom), geom; kw...)
 # There is no trait and this is an AbstractArray - so just iterate over it calling _apply on the contents
 @inline function _apply(f::F, target, ::Nothing, A::AbstractArray; threaded, kw...) where F
-    # For an Array there is nothing else to do but map `_apply` over all values
-    # _maptasks may run this level threaded if `threaded==true`,
-    # but deeper `_apply` called in the closure will not be threaded
-    apply_to_array(i) = _apply(f, target, A[i]; threaded=False(), kw...)
-    _maptasks(apply_to_array, eachindex(A), threaded)
+    applicator = ApplyToArray(f, target, A; kw...)
+    _maptasks(applicator, eachindex(A), threaded)
 end
 # There is no trait and this is not an AbstractArray.
 # Try to call _apply over it. We can't use threading
 # as we don't know if we can can index into it. So just `map`.
 @inline function _apply(f::F, target, ::Nothing, iterable::IterableType; threaded, kw...) where {F, IterableType}
     # Try the Tables.jl interface first
     if Tables.istable(iterable)
-    _apply_table(f, target, iterable; threaded, kw...)
+        _apply_table(f, target, iterable; threaded, kw...)
     else # this is probably some form of iterable...
         if threaded isa True
             # `collect` first so we can use threads
@@ -239,9 +236,8 @@ end
 ) where F
 
     # Run _apply on all `features` in the feature collection, possibly threaded
-    apply_to_feature(i) =
-        _apply(f, target, GI.getfeature(fc, i); crs, calc_extent, threaded=False())::GI.Feature
-    features = _maptasks(apply_to_feature, 1:GI.nfeature(fc), threaded)
+    applicator = ApplyToFeatures(f, target, fc; crs, calc_extent)
+    features = _maptasks(applicator, 1:GI.nfeature(fc), threaded)
     if calc_extent isa True
         # Calculate the extent of the features
         extent = mapreduce(GI.extent, Extents.union, features)
@@ -278,21 +274,16 @@ end
     # Map `_apply` over all sub geometries of `geom`
     # to create a new vector of geometries
     # TODO handle zero length
-    apply_to_geom(i) = _apply(f, target, GI.getgeom(geom, i); crs, calc_extent, threaded=False())
-    geoms = _maptasks(apply_to_geom, 1:GI.ngeom(geom), threaded)
+    applicator = ApplyToGeom(f, target, geom; crs, calc_extent)
+    geoms = _maptasks(applicator, 1:GI.ngeom(geom), threaded)
     return _apply_inner(geom, geoms, crs, calc_extent)
 end
 @inline function _apply(f::F, target::TraitTarget{<:PointTrait}, trait::GI.PolygonTrait, geom;
     crs=GI.crs(geom), calc_extent=False(), threaded
 )::(GI.geointerface_geomtype(trait)) where F
     # We need to force rebuilding a LinearRing not a LineString
-    geoms = _maptasks(1:GI.ngeom(geom), threaded) do i
-        lr = GI.getgeom(geom, i)
-        points = map(GI.getgeom(lr)) do p
-            _apply(f, target, p; crs, calc_extent, threaded=False())
-        end
-        _linearring(_apply_inner(lr, points, crs, calc_extent))
-    end
+    applicator = ApplyPointsToPolygon(f, target, geom; crs, calc_extent)
+    geoms = _maptasks(applicator, 1:GI.ngeom(geom), threaded)
     return _apply_inner(geom, geoms, crs, calc_extent)
 end
 function _apply_inner(geom, geoms, crs, calc_extent::True)
@@ -360,6 +351,24 @@ end
         end
     end
 
+    # Finally we join the results into a new vector
+    return mapreduce(fetch, vcat, tasks)
+end
+@inline function _maptasks(a::Applicator{<:TaskFunctors}, taskrange, threaded::True)::Vector
+    ntasks = length(taskrange)
+    chunk_size = max(1, cld(ntasks, (length(a.f.functors))))
+    # partition the range into chunks
+    task_chunks = Iterators.partition(taskrange, chunk_size)
+    # Map over the chunks
+    tasks = map(task_chunks, view(a.f.functors, 1:length(task_chunks))) do chunk, ft
+        f = rebuild(a, ft)
+        # Spawn a task to process this chunk
+        StableTasks.@spawn begin
+            # Where we map `f` over the chunk indices
+            map(f, chunk)
+        end
+    end
+
     # Finally we join the results into a new vector
     return mapreduce(fetch, vcat, tasks)
 end

GeometryOpsCore/src/types/applicators.jl

@@ -0,0 +1,108 @@
+#=
+
+```@meta
+CollapsedDocStrings = true
+```
+# Applicators
+
+Applicators are functor structs that wrap a function and a target object.  They are used to
+implement the `apply` interface, and are also used to dispatch to the correct method of `_maptasks`.
+
+The basic functor struct is `ApplyToGeom`, which just applies a function to a geometry.
+
+```@docs; canonical=false
+Applicator
+ApplyToGeom
+ApplyToArray
+ApplyToFeatures
+```
+=#
+
+"""
+    abstract type Applicator{F,T}
+
+An abstract type for applicators that apply a function to a target object.
+
+The type parameter `F` is the type of the function to apply, and `T` is the type of the target object.
+
+A common dispatch pattern is to dispatch on `F` which may also be e.g. a ThreadFunctor.
+
+## Interface
+All applicators must be callable by an index integer, and define the following methods:
+
+- `rebuild(a::Applicator, f)` - swap out the function and return a new applicator.
+
+The calling convention is `my_applicator(i::Int)`, so applicators must define this method.
+"""
+abstract type Applicator{F,T} end
+
+struct ApplyToPoint{Z,M,F} <: Applicator{F,Nothing} 
+    f::F
+end
+ApplyToPoint{Z,M}(f::F) where {Z,M,F} = ApplyToPoint{Z,M,F}(f)
+ApplyToPoint{Z}(f::F) where {Z,F} = ApplyToPoint{Z,false,F}(f)
+# Default function is just `tuple`
+(a::Type{<:ApplyToPoint})() = a(tuple)
+rebuild(a::ApplyToPoint{Z, M}, f::F) where {Z, M, F} = ApplyToPoint{Z, M, F}(f)
+
+# Currently we ignore M by default
+const WithXY = ApplyToPoint{false}
+const WithXYZ = ApplyToPoint{true}
+# But these could be used to require M
+const WithXYM = ApplyToPoint{false,true}
+const WithXYZM = ApplyToPoint{true,true}
+
+(t::WithXY)(p) = t.f(GI.x(p), GI.y(p))
+(t::WithXYZ)(p) = t.f(GI.x(p), GI.y(p), GI.z(p))
+(t::WithXYZM)(p) = t.f(GI.x(p), GI.y(p), GI.m(p))
+(t::WithXYM)(p) = t.f(GI.x(p), GI.y(p), GI.z(p), GI.m(p))
+
+# ***
+
+for T in (:ApplyToGeom, :ApplyToArray, :ApplyToFeatures, :ApplyPointsToPolygon)
+    @eval begin
+        struct $T{F,T,O,K} <: Applicator{F,T}
+            f::F
+            target::T
+            obj::O
+            kw::K
+        end
+        $T(f, target, obj; kw...) = $T(f, target, obj, kw)
+        # rebuild lets us swap out the function, such as with TaskFunctors
+        rebuild(a::$T, f) = $T(f, a.target, a.obj, a.kw) 
+    end
+end
+
+# Functor definitions
+# _maptasks may run this level threaded if `threaded==true`
+# but deeper `_apply` calls will not be threaded
+# For an Array there is nothing to do but map `_apply` over all values
+(a::ApplyToArray)(i::Int) = _apply(a.f, a.target, a.obj[i]; a.kw..., threaded=False())
+# For a FeatureCollection or Geometry we need getfeature or getgeom calls
+(a::ApplyToFeatures)(i::Int) = _apply(a.f, a.target, GI.getfeature(a.obj, i); a.kw..., threaded=False())
+(a::ApplyToGeom)(i::Int) = _apply(a.f, a.target, GI.getgeom(a.obj, i); a.kw..., threaded=False())
+function (a::ApplyPointsToPolygon)(i::Int)
+    lr = GI.getgeom(a.obj, i)
+    points = map(GI.getgeom(lr)) do p
+        _apply(a.f, a.target, p; a.kw..., threaded=False())
+    end
+    _linearring(_apply_inner(lr, points, a.kw[:crs], a.kw[:calc_extent]))
+end
+
+@doc """
+    ApplyToArray(f, target, arr, kw)
+
+Create an [`Applicator`](@ref) that applies a function to all elements of `arr`.
+""" ApplyToArray
+
+@doc """
+    ApplyToGeom(f, target, geom, kw)
+
+Create an [`Applicator`](@ref) that applies a function to all sub-geometries of `geom`.
+""" ApplyToGeom
+
+@doc """
+    ApplyToFeatures(f, target, fc, kw)
+
+Create an [`Applicator`](@ref) that applies a function to all features of `fc`.
+""" ApplyToFeatures

GeometryOpsCore/src/types/booltypes.jl

@@ -51,3 +51,7 @@ function booltype end
 
 @inline booltype(x::Bool)::BoolsAsTypes = x ? True() : False()
 @inline booltype(x::BoolsAsTypes)::BoolsAsTypes = x
+
+@inline istrue(x::True) = true
+@inline istrue(x::False) = false
+@inline istrue(x::Bool) = x

GeometryOpsCore/src/types/tasks.jl

@@ -0,0 +1,35 @@
+#=
+# TaskFunctors
+
+Many functions a user might run through `apply` call into C.  But what do you do
+when the C function you are calling is not threadsafe, or uses a reentrant/context
+based API?
+
+The answer is to wrap the function in TaskFunctors, which will use an approximation of
+task local storage to ensure that each task calls its own copy of the function or C object
+you are invoking.
+
+The primary application of this is Proj `reproject`, and Proj has a reentrant API based on context
+objects.  So, to make this work, we need to create a new context object per task.  We do this and 
+pass the vector of Proj transformation objects to `TaskFunctors`, which `apply` and `_maptasks`
+dispatch on to behave correctly in this circumstance.
+
+```@docs; canonical=false
+TaskFunctors
+```
+=#
+"""
+    TaskFunctors(functors, tasks_per_thread)
+
+A struct to hold the functors and tasks_per_thread,
+for internal use with `_maptasks` where functions have state
+that cannot be shared accross threads, such as `Proj.Transformation`.
+
+`functors` must be an array or tuple of functors, one per thread,
+and `tasks_per_thread` must be an integer.  This also allows you to control
+the number of tasks per thread that `_maptasks` launches, useful for tuning 
+performance if you like.
+"""
+struct TaskFunctors{F} 
+    functors::F
+end

Project.toml

@@ -8,6 +8,7 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
+GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 GeometryOpsCore = "05efe853-fabf-41c8-927e-7063c8b9f013"
@@ -32,6 +33,7 @@ DataAPI = "1"
 DelaunayTriangulation = "1.0.4"
 ExactPredicates = "2.2.8"
 FlexiJoins = "0.1.30"
+GeoFormatTypes = "0.4"
 GeoInterface = "1.2"
 GeometryBasics = "0.4.7, 0.5"
 GeometryOpsCore = "=0.1.3"

ext/GeometryOpsProjExt/reproject.jl

@@ -1,6 +1,12 @@
-import GeometryOps: GI, GeoInterface, reproject, apply, transform, _is3d, True, False
+import GeometryOps: GI, GeoInterface, reproject, apply, transform, _is3d, istrue,
+    True, False, TaskFunctors, WithXY, WithXYZ
+import GeoFormatTypes
 import Proj
 
+# TODO:
+# - respect `time`
+# - respect measured values
+
 function reproject(geom;
     source_crs=nothing, target_crs=nothing, transform=nothing, kw...
 )
@@ -22,29 +28,53 @@ function reproject(geom;
         reproject(geom, transform; kw...)
     end
 end
-function reproject(geom, source_crs, target_crs;
-    time=Inf,
-    always_xy=true,
-    transform=nothing,
+function reproject(geom, source_crs, target_crs; always_xy=true, kw...)
+    transform = Proj.Transformation(convert(String, source_crs), convert(String, target_crs); always_xy)
+    return reproject(geom, transform; target_crs, kw...)
+end
+function reproject(
+    geom, source_crs::CRSType, target_crs::CRSType; always_xy=true, kw...
+) where CRSType <: Union{GeoFormatTypes.GeoFormat, Proj.CRS}
+    transform = Proj.Transformation(source_crs, target_crs; always_xy)
+    return reproject(geom, transform; target_crs, kw...)
+end
+function reproject(geom, transform::Proj.Transformation; 
+    context=C_NULL, 
+    target_crs=nothing, 
+    time=Inf, 
+    threaded=False(), 
     kw...
 )
-    transform = if isnothing(transform) 
-        s = source_crs isa Proj.CRS ? source_crs : convert(String, source_crs)
-        t = target_crs isa Proj.CRS ? target_crs : convert(String, target_crs)
-        Proj.Transformation(s, t; always_xy)
-    else
-        transform
+    if isnothing(target_crs)
+        target_crs = GeoFormatTypes.ESRIWellKnownText(Proj.CRS(Proj.proj_get_target_crs(transform.pj)))
     end
-    reproject(geom, transform; time, target_crs, kw...)
-end
-function reproject(geom, transform::Proj.Transformation; time=Inf, target_crs=nothing, kw...)
-    if _is3d(geom)
-        return apply(GI.PointTrait(), geom; crs=target_crs, kw...) do p
-            transform(GI.x(p), GI.y(p), GI.z(p))
+    kw1 = (; crs=target_crs, threaded, kw...)
+    if istrue(threaded)
+        tasks_per_thread = 2
+        ntasks = Threads.nthreads() * tasks_per_thread
+        # Construct one context per planned task
+        contexts = [Proj.proj_context_clone(context) for _ in 1:ntasks]
+        # Clone the transformation for each context
+        proj_transforms = [Proj.Transformation(Proj.proj_clone(transform.pj)) for context in contexts]
+        # Assign the context to the transformation
+        Proj.proj_assign_context.(getproperty.(proj_transforms, :pj), contexts)
+
+        results = if _is3d(geom)
+            functors = TaskFunctors(WithXYZ.(proj_transforms))
+            apply(functors, GI.PointTrait(), geom; kw1...)
+        else
+            functors = TaskFunctors(WithXY.(proj_transforms))
+            apply(functors, GI.PointTrait(), geom; kw1...)
         end
+        # Destroy the temporary threading contexts that we created
+        Proj.proj_destroy.(contexts)
+        # Return the results
+        return results
     else
-        return apply(GI.PointTrait(), geom; crs=target_crs, kw...) do p
-            transform(GI.x(p), GI.y(p))
+        if _is3d(geom)
+            return apply(WithXYZ(transform), GI.PointTrait(), geom; kw1...)
+        else
+            return apply(WithXY(transform), GI.PointTrait(), geom; kw1...)
         end
-    end
-end
+    end    
+end

src/GeometryOps.jl

@@ -7,7 +7,9 @@ import GeometryOpsCore:
                 TraitTarget,
                 Manifold, Planar, Spherical, Geodesic, AutoManifold, WrongManifoldException,
                 Algorithm, AutoAlgorithm, ManifoldIndependentAlgorithm, SingleManifoldAlgorithm, NoAlgorithm,
-                BoolsAsTypes, True, False, booltype,
+                BoolsAsTypes, True, False, booltype, istrue,
+                TaskFunctors, 
+                WithXY, WithXYZ, WithXYM, WithXYZM,
                 apply, applyreduce, 
                 flatten, reconstruct, rebuild, unwrap, _linearring,
                 APPLY_KEYWORDS, THREADED_KEYWORD, CRS_KEYWORD, CALC_EXTENT_KEYWORD

src/transformations/forcedims.jl

@@ -13,11 +13,7 @@ export forcexy, forcexyz
 
 Force the geometry to be 2D.  Works on any geometry, vector of geometries, feature collection, or table!
 """
-function forcexy(geom)
-    return apply(GI.PointTrait(), geom) do point
-        (GI.x(point), GI.y(point))
-    end
-end
+forcexy(geom) = apply(WithXY(), GI.PointTrait(), geom)
 
 """
     forcexyz(geom, z = 0)