speed up large polygon intersections using STRtrees and many other things

.github/workflows/CI.yml

@@ -18,7 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9'
+          - '1.10'
           - '1'
           - 'nightly'
         os:
@@ -33,8 +33,8 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v1
-      - name: Dev GeometryOpsCore`
-        run: julia --project=. -e 'using Pkg; Pkg.develop(; path = joinpath(".", "GeometryOpsCore"))'
+      - name: Dev GeometryOpsCore and add other packages
+        run: julia --project=. -e 'using Pkg; Pkg.develop(; path = joinpath(".", "GeometryOpsCore"));'
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
@@ -60,13 +60,28 @@ jobs:
         with:
           version: '1'
       - name: Build and add versions
-        run: julia --project=docs -e 'using Pkg; Pkg.develop([PackageSpec(path = "."), PackageSpec(path = joinpath(".", "GeometryOpsCore"))]); Pkg.add([PackageSpec(name = "GeoMakie", rev = "master"), PackageSpec(name="GeoInterface", rev="bugfix_vars")])'
+        run: julia --project=docs -e 'using Pkg; Pkg.add([PackageSpec(name = "GeoMakie", rev = "master")])'
       - uses: julia-actions/julia-docdeploy@v1
         with:
           install-package: false
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
+  doctests:
+    name: Doctests
+    runs-on: ubuntu-latest
+    permissions:
+      contents: write
+      statuses: write
+      actions: write
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/cache@v1
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: '1'
+      - name: Build and add versions
+        run: julia --project=docs -e 'using Pkg; Pkg.add([PackageSpec(name = "GeoMakie", rev = "master")])'
       - run: |
           julia --project=docs -e '
             using Documenter: DocMeta, doctest

GeometryOpsCore/src/GeometryOpsCore.jl

@@ -19,6 +19,7 @@ using DataAPI
 import StableTasks
 
 include("keyword_docs.jl")
+include("constants.jl")
 
 include("types/manifold.jl")
 include("types/algorithm.jl")

GeometryOpsCore/src/apply.jl

@@ -126,8 +126,8 @@ end
 @inline function apply(
     f::F, target, geom; calc_extent=false, threaded=false, kw...
 ) where F
-    threaded = _booltype(threaded)
-    calc_extent = _booltype(calc_extent)
+    threaded = booltype(threaded)
+    calc_extent = booltype(calc_extent)
     _apply(f, TraitTarget(target), geom; threaded, calc_extent, kw...)
 end
 
@@ -139,7 +139,7 @@ end
     # 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...)
+    apply_to_array(i) = _apply(f, target, A[i]; threaded=False(), kw...)
     _maptasks(apply_to_array, eachindex(A), threaded)
 end
 # There is no trait and this is not an AbstractArray.
@@ -150,7 +150,7 @@ end
     if Tables.istable(iterable)
     _apply_table(f, target, iterable; threaded, kw...)
     else # this is probably some form of iterable...
-        if threaded isa _True
+        if threaded isa True
             # `collect` first so we can use threads
             _apply(f, target, collect(iterable); threaded, kw...)
         else
@@ -226,14 +226,14 @@ end
 # Rewrap all FeatureCollectionTrait feature collections as GI.FeatureCollection
 # Maybe use threads to call _apply on component features
 @inline function _apply(f::F, target, ::GI.FeatureCollectionTrait, fc;
-    crs=GI.crs(fc), calc_extent=_False(), threaded
+    crs=GI.crs(fc), calc_extent=False(), threaded
 ) 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
+        _apply(f, target, GI.getfeature(fc, i); crs, calc_extent, threaded=False())::GI.Feature
     features = _maptasks(apply_to_feature, 1:GI.nfeature(fc), threaded)
-    if calc_extent isa _True
+    if calc_extent isa True
         # Calculate the extent of the features
         extent = mapreduce(GI.extent, Extents.union, features)
         # Return a FeatureCollection with features, crs and calculated extent
@@ -245,13 +245,13 @@ end
 end
 # Rewrap all FeatureTrait features as GI.Feature, keeping the properties
 @inline function _apply(f::F, target, ::GI.FeatureTrait, feature;
-    crs=GI.crs(feature), calc_extent=_False(), threaded
+    crs=GI.crs(feature), calc_extent=False(), threaded
 ) where F
     # Run _apply on the contained geometry
     geometry = _apply(f, target, GI.geometry(feature); crs, calc_extent, threaded)
     # Get the feature properties
     properties = GI.properties(feature)
-    if calc_extent isa _True
+    if calc_extent isa True
         # Calculate the extent of the geometry
         extent = GI.extent(geometry)
         # Return a new Feature with the new geometry and calculated extent, but the original properties and crs
@@ -264,36 +264,36 @@ end
 # Reconstruct nested geometries,
 # maybe using threads to call _apply on component geoms
 @inline function _apply(f::F, target, trait, geom;
-    crs=GI.crs(geom), calc_extent=_False(), threaded
+    crs=GI.crs(geom), calc_extent=False(), threaded
 )::(GI.geointerface_geomtype(trait)) where F
     # 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())
+    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)
     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
+    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())
+            _apply(f, target, p; crs, calc_extent, threaded=False())
         end
         _linearring(_apply_inner(lr, points, crs, calc_extent))
     end
     return _apply_inner(geom, geoms, crs, calc_extent)
 end
-function _apply_inner(geom, geoms, crs, calc_extent::_True)
+function _apply_inner(geom, geoms, crs, calc_extent::True)
     # Calculate the extent of the sub geometries
     extent = mapreduce(GI.extent, Extents.union, geoms)
     # Return a new geometry of the same trait as `geom`,
     # holding the new `geoms` with `crs` and calculated extent
     return rebuild(geom, geoms; crs, extent)
 end
-function _apply_inner(geom, geoms, crs, calc_extent::_False)
+function _apply_inner(geom, geoms, crs, calc_extent::False)
     # Return a new geometry of the same trait as `geom`, holding the new `geoms` with `crs`
     return rebuild(geom, geoms; crs)
 end
@@ -319,41 +319,22 @@ end
 
 using Base.Threads: nthreads, @threads, @spawn
 
-# Threading utility, modified Mason Protters threading PSA
-# run `f` over ntasks, where f receives an AbstractArray/range
-# of linear indices
-@inline function _maptasks(f::F, taskrange, threaded::_True)::Vector where F
-    ntasks = length(taskrange)
-    # Customize this as needed.
-    # More tasks have more overhead, but better load balancing
-    tasks_per_thread = 2
-    chunk_size = max(1, ntasks ÷ (tasks_per_thread * nthreads()))
-    # partition the range into chunks
-    task_chunks = Iterators.partition(taskrange, chunk_size)
-    # Map over the chunks
-    tasks = map(task_chunks) do chunk
-        # Spawn a task to process this chunk
-        @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
 #=
 Here we used to use the compiler directive `@assume_effects :foldable` to force the compiler
 to lookup through the closure. This alone makes e.g. `flip` 2.5x faster!
 
 But it caused inference to fail, so we've removed it.  No effect on runtime so far as we can tell, 
 at least in Julia 1.11.
 =#
-@inline function _maptasks(f::F, taskrange, threaded::_False)::Vector where F
+@inline function _maptasks(f::F, taskrange, threaded::False)::Vector where F
     map(f, taskrange)
 end
 
-@inline function _maptasks(f::F, taskrange, threaded::_TrueButStable)::Vector where F
+
+# Threading utility, modified Mason Protters threading PSA
+# run `f` over ntasks, where f receives an AbstractArray/range
+# of linear indices
+@inline function _maptasks(f::F, taskrange, threaded::True)::Vector where F
     ntasks = length(taskrange)
     # Customize this as needed.
     # More tasks have more overhead, but better load balancing

GeometryOpsCore/src/applyreduce.jl

@@ -34,24 +34,24 @@ feature collections and nested geometries.
 @inline function applyreduce(
     f::F, op::O, target, geom; threaded=false, init=nothing
 ) where {F, O}
-    threaded = _booltype(threaded)
+    threaded = booltype(threaded)
     _applyreduce(f, op, TraitTarget(target), geom; threaded, init)
 end
 
 @inline _applyreduce(f::F, op::O, target, geom; threaded, init) where {F, O} =
     _applyreduce(f, op, target, GI.trait(geom), geom; threaded, init)
 # Maybe use threads reducing over arrays
 @inline function _applyreduce(f::F, op::O, target, ::Nothing, A::AbstractArray; threaded, init) where {F, O}
-    applyreduce_array(i) = _applyreduce(f, op, target, A[i]; threaded=_False(), init)
+    applyreduce_array(i) = _applyreduce(f, op, target, A[i]; threaded=False(), init)
     _mapreducetasks(applyreduce_array, op, eachindex(A), threaded; init)
 end
 # Try to applyreduce over iterables
 @inline function _applyreduce(f::F, op::O, target, ::Nothing, iterable::IterableType; threaded, init) where {F, O, IterableType}
     if Tables.istable(iterable)
         _applyreduce_table(f, op, target, iterable; threaded, init)
     else
-        applyreduce_iterable(i) = _applyreduce(f, op, target, i; threaded=_False(), init)
-        if threaded isa _True # Try to `collect` and reduce over the vector with threads
+        applyreduce_iterable(i) = _applyreduce(f, op, target, i; threaded=False(), init)
+        if threaded isa True # Try to `collect` and reduce over the vector with threads
             _applyreduce(f, op, target, collect(iterable); threaded, init)
         else
             # Try to `mapreduce` the iterable as-is
@@ -77,23 +77,23 @@ function _applyreduce_table(f::F, op::O, target::GI.FeatureTrait, iterable::Iter
 end
 # Maybe use threads reducing over features of feature collections
 @inline function _applyreduce(f::F, op::O, target, ::GI.FeatureCollectionTrait, fc; threaded, init) where {F, O}
-    applyreduce_fc(i) = _applyreduce(f, op, target, GI.getfeature(fc, i); threaded=_False(), init)
+    applyreduce_fc(i) = _applyreduce(f, op, target, GI.getfeature(fc, i); threaded=False(), init)
     _mapreducetasks(applyreduce_fc, op, 1:GI.nfeature(fc), threaded; init)
 end
 # Features just applyreduce to their geometry
 @inline _applyreduce(f::F, op::O, target, ::GI.FeatureTrait, feature; threaded, init) where {F, O} =
     _applyreduce(f, op, target, GI.geometry(feature); threaded, init)
 # Maybe use threads over components of nested geometries
 @inline function _applyreduce(f::F, op::O, target, trait, geom; threaded, init) where {F, O}
-    applyreduce_geom(i) = _applyreduce(f, op, target, GI.getgeom(geom, i); threaded=_False(), init)
+    applyreduce_geom(i) = _applyreduce(f, op, target, GI.getgeom(geom, i); threaded=False(), init)
     _mapreducetasks(applyreduce_geom, op, 1:GI.ngeom(geom), threaded; init)
 end
 # Don't thread over points it won't pay off
 @inline function _applyreduce(
     f::F, op::O, target, trait::Union{GI.LinearRing,GI.LineString,GI.MultiPoint}, geom;
     threaded, init
 ) where {F, O}
-    _applyreduce(f, op, target, GI.getgeom(geom); threaded=_False(), init)
+    _applyreduce(f, op, target, GI.getgeom(geom); threaded=False(), init)
 end
 # Apply f to the target
 @inline function _applyreduce(f::F, op::O, ::TraitTarget{Target}, ::Trait, x; kw...) where {F,O,Target,Trait<:Target} 
@@ -124,29 +124,7 @@ import Base.Threads: nthreads, @threads, @spawn
 # 
 # If you absolutely need a single chunk, then `threaded = false` will always decompose
 # to straight `mapreduce` without grouping.
-
-@inline function _mapreducetasks(f::F, op, taskrange, threaded::_True; init) where F
-    ntasks = length(taskrange)
-    # Customize this as needed.
-    # More tasks have more overhead, but better load balancing
-    tasks_per_thread = 2
-    chunk_size = max(1, ntasks ÷ (tasks_per_thread * nthreads()))
-    # partition the range into chunks
-    task_chunks = Iterators.partition(taskrange, chunk_size)
-    # Map over the chunks
-    tasks = map(task_chunks) do chunk
-        # Spawn a task to process this chunk
-        Threads.@spawn begin
-            # Where we map `f` over the chunk indices
-            mapreduce(f, op, chunk; init)
-        end
-    end
-
-    # Finally we join the results into a new vector
-    return mapreduce(fetch, op, tasks; init)
-end
-
-@inline function _mapreducetasks(f::F, op, taskrange, threaded::_TrueButStable; init) where F
+@inline function _mapreducetasks(f::F, op, taskrange, threaded::True; init) where F
     ntasks = length(taskrange)
     # Customize this as needed.
     # More tasks have more overhead, but better load balancing
@@ -167,6 +145,6 @@ end
     return mapreduce(fetch, op, tasks; init)
 end
 
-function _mapreducetasks(f::F, op, taskrange, threaded::_False; init) where F
+function _mapreducetasks(f::F, op, taskrange, threaded::False; init) where F
     mapreduce(f, op, taskrange; init)
 end

GeometryOpsCore/src/constants.jl

@@ -0,0 +1,8 @@
+"The semi-major axis of the WGS84 ellipsoid"
+const WGS84_EARTH_SEMI_MAJOR_RADIUS = 6378137.0
+
+"The inverse flattening of the WGS84 ellipsoid"
+const WGS84_EARTH_INV_FLATTENING = 298.257223563
+
+"The mean radius of the WGS84 ellipsoid, used for spherical manifold default"
+const WGS84_EARTH_MEAN_RADIUS = 6371008.8

GeometryOpsCore/src/types/algorithm.jl

@@ -14,8 +14,19 @@ Algorithms are:
 
 =#
 
+export Algorithm, AutoAlgorithm, ManifoldIndependentAlgorithm, SingleManifoldAlgorithm, NoAlgorithm
+
 abstract type Algorithm{M <: Manifold} end
 
+struct AutoAlgorithm{T, M <: Manifold} <: Algorithm{M} 
+    manifold::M
+    x::T
+end
+
+AutoAlgorithm(m::Manifold; kwargs...) = AutoAlgorithm(m, kwargs)
+AutoAlgorithm(; kwargs...) = AutoAlgorithm(AutoManifold(), kwargs)
+
+
 abstract type ManifoldIndependentAlgorithm{M <: Manifold} <: Algorithm{M} end
 
 abstract type SingleManifoldAlgorithm{M <: Manifold} <: Algorithm{M} end