[Core] Add an algorithm interface + StableTasks (#269)

This was factored out of the "dev branch" #259 and contains the subset of changes that apply to GeometryOpsCore, for easier review.

Child PRs: #271 (TGGeometry) -> #275 (AdaptivePredicates) -> #273 (clipping algorithm type) -> #274 (trees)

- Use [StableTasks.jl](https://github.com/JuliaFolds2/StableTasks.jl) in apply and applyreduce - its type-stable tasks save us some allocations!
- Remove `Base.@assume_effects` on the low level functions, which caused issues on Julia v1.11 and was probably incorrect anyway
- Add an algorithm interface with an abstract supertype `Algorithm{M <: Manifold}`, as discussed in #247.  Also adds an abstract Operator supertype and some discussion in code comments, but no implementation or interface surface there yet.
- Split out `types.jl` into a directory `types` with a bunch of files in it, for ease of readability / docs / use.
- (out of context change): refactor CI a bit for cleanliness.


TODOs for later (not this PR):
- [ ] Add a `format` method that takes in an incompletely specified algorithm and some geometry as input, and returns a completely specified algorithm.  What does this mean?  Imagine I call `GO.intersection(FosterHormannClipping(), geom1, geom2)`.  That `FosterHormannClipping()` should get expanded to `FosterHormannClipping(AutoAlgorithm(), AutoAccelerator())`.  Then, `format` will take `format(alg, args...)` and:
  - get the `crstrait` of the two geometries, scan for incompatibilities, assign the correct manifold to the algorithm (maybe warn or emit debug info)
  - if no geometries available, get the manifold via `best_manifold(::Algorithm)`.
  - maybe inflate the accelerator by checking `npoint` and later preparations to see what's most efficient, maybe not - depends on what we want!

Author: asinghvi17

.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

.gitignore

@@ -3,4 +3,6 @@
 /docs/build/
 /docs/src/source/
 .vscode/
-.DS_Store
+.DS_Store
+
+benchmarks/Manifest.toml

GeometryOpsCore/Project.toml

@@ -1,15 +1,17 @@
 name = "GeometryOpsCore"
 uuid = "05efe853-fabf-41c8-927e-7063c8b9f013"
 authors = ["Anshul Singhvi <[email protected]>", "Rafael Schouten <[email protected]>", "Skylar Gering <[email protected]>", "and contributors"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
+StableTasks = "91464d47-22a1-43fe-8b7f-2d57ee82463f"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
-julia = "1.9"
 DataAPI = "1"
 GeoInterface = "1.2"
+StableTasks = "0.1.5"
 Tables = "1"
+julia = "1.9"

GeometryOpsCore/README.md

@@ -2,7 +2,13 @@
 
 This is a "core" package for [GeometryOps.jl](https://github.com/JuliaGeo/GeometryOps.jl), that defines some basic primitive functions and types for GeometryOps.
 
-Generally, you would depend on this to use either the GeometryOps types (like `Linear`, `Spherical`, etc) or the primitive functions like `apply`, `applyreduce`, `flatten`, etc. 
+It defines, all in all:
+- Manifolds and the manifold interface
+- The Algorithm type and the algorithm interface
+- Low level functions like apply, applyreduce, flatten, etc.
+- Common methods that should work across all geometries!
+
+Generally, you would depend on this to use either the GeometryOps types (like `Planar`, `Spherical`, etc) or the primitive functions like `apply`, `applyreduce`, `flatten`, etc. 
 All of these are also accessible from GeometryOps, so it's preferable that you use GeometryOps directly.
 
 Tests are in the main GeometryOps tests, we don't have separate tests for GeometryOpsCore since it's in a monorepo structure.

GeometryOpsCore/src/GeometryOpsCore.jl

@@ -6,7 +6,7 @@ import GeoInterface
 import GeoInterface as GI
 import GeoInterface: Extents
 
-# Import all names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`.
+# Import all exported names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`.
 for name in names(GeoInterface)
     @eval using GeoInterface: $name
 end
@@ -16,9 +16,17 @@ end
 
 using Tables
 using DataAPI
+import StableTasks
 
 include("keyword_docs.jl")
-include("types.jl")
+include("constants.jl")
+
+include("types/manifold.jl")
+include("types/algorithm.jl")
+include("types/operation.jl")
+include("types/exceptions.jl")
+include("types/booltypes.jl")
+include("types/traittarget.jl")
 
 include("apply.jl")
 include("applyreduce.jl")

GeometryOpsCore/src/apply.jl

@@ -30,18 +30,11 @@ Functions like [`flip`](@ref), [`reproject`](@ref), [`transform`](@ref), even [`
 using the `apply` framework.  Similarly, [`centroid`](@ref), [`area`](@ref) and [`distance`](@ref) have been implemented using the 
 [`applyreduce`](@ref) framework.
 
-## Docstrings
-
-### Functions
-
 ```@docs; collapse=true, canonical=false
 apply
-applyreduce
 ```
 
-=#
 
-#=
 ## What is `apply`?
 
 `apply` applies some function to every geometry matching the `Target`
@@ -69,7 +62,7 @@ Be careful making a union across "levels" of nesting, e.g.
 `Union{FeatureTrait,PolygonTrait}`, as `_apply` will just never reach
 `PolygonTrait` when all the polygons are wrapped in a `FeatureTrait` object.
 
-## Embedding:
+### Embedding
 
 `extent` and `crs` can be embedded in all geometries, features, and
 feature collections as part of `apply`. Geometries deeper than `Target`
@@ -78,14 +71,30 @@ will of course not have new `extent` or `crs` embedded.
 - `calc_extent` signals to recalculate an `Extent` and embed it.
 - `crs` will be embedded as-is
 
-## Threading
+### Threading
 
 Threading is used at the outermost level possible - over
 an array, feature collection, or e.g. a MultiPolygonTrait where
 each `PolygonTrait` sub-geometry may be calculated on a different thread.
 
 Currently, threading defaults to `false` for all objects, but can be turned on
 by passing the keyword argument `threaded=true` to `apply`.
+
+Threading uses [StableTasks.jl](https://github.com/JuliaFolds2/StableTasks.jl) to provide
+type-stable tasks (base Julia `Threads.@spawn` is not type stable).  This is completely cost-free
+and saves some allocations when running multithreaded. 
+
+The current strategy is to launch 2 tasks for each CPU thread, to provide load balancing.  We
+assume Julia will manage these tasks efficiently, and we don't want to run too many tasks
+since each task does have some overhead when it's created.  This may need revisiting in the future,
+but it's a pretty easy heuristic to use.
+
+## Implementation
+
+Literate.jl source code is below.
+
+***
+
 =#
 
 """
@@ -319,6 +328,18 @@ end
 
 using Base.Threads: nthreads, @threads, @spawn
 
+#=
+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
+    map(f, taskrange)
+end
+
+
 # Threading utility, modified Mason Protters threading PSA
 # run `f` over ntasks, where f receives an AbstractArray/range
 # of linear indices
@@ -333,22 +354,12 @@ using Base.Threads: nthreads, @threads, @spawn
     # Map over the chunks
     tasks = map(task_chunks) do chunk
         # Spawn a task to process this chunk
-        @spawn begin
+        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
-#=
-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
-    map(f, taskrange)
 end

GeometryOpsCore/src/applyreduce.jl

@@ -18,6 +18,36 @@ and perform some operation on it.
 
 [`centroid`](@ref), [`area`](@ref) and [`distance`](@ref) have been implemented using the 
 [`applyreduce`](@ref) framework.
+
+```@docs
+applyreduce
+```
+
+
+### Threading
+
+Threading is used at the outermost level possible - over
+an array, feature collection, or e.g. a MultiPolygonTrait where
+each `PolygonTrait` sub-geometry may be calculated on a different thread.
+
+Currently, threading defaults to `false` for all objects, but can be turned on
+by passing the keyword argument `threaded=true` to `apply`.
+
+Threading uses [StableTasks.jl](https://github.com/JuliaFolds2/StableTasks.jl) to provide
+type-stable tasks (base Julia `Threads.@spawn` is not type stable).  This is completely cost-free
+and saves some allocations when running multithreaded. 
+
+The current strategy is to launch 2 tasks for each CPU thread, to provide load balancing.  We
+assume Julia will manage these tasks efficiently, and we don't want to run too many tasks
+since each task does have some overhead when it's created.  This may need revisiting in the future,
+but it's a pretty easy heuristic to use.
+
+## Implementation
+
+Literate.jl source code is below.
+
+***
+
 =#
 
 """
@@ -135,7 +165,7 @@ import Base.Threads: nthreads, @threads, @spawn
     # Map over the chunks
     tasks = map(task_chunks) do chunk
         # Spawn a task to process this chunk
-        @spawn begin
+        StableTasks.@spawn begin
             # Where we map `f` over the chunk indices
             mapreduce(f, op, chunk; init)
         end
@@ -144,6 +174,7 @@ import Base.Threads: nthreads, @threads, @spawn
     # Finally we join the results into a new vector
     return mapreduce(fetch, op, tasks; init)
 end
-Base.@assume_effects :foldable 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