Rewrite for version 0.2

.github/workflows/CI.yml

@@ -18,39 +18,21 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.8'
-          - '1.9'
-          - 'nightly'
+          - 'pre'
+          - 'min'
         os:
           - ubuntu-latest
-          - windows-latest
-          - macos-latest
         arch:
-          - x64
+          - 'default'
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v7
+      - uses: julia-actions/setup-julia@v3
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v1
-      - name: Develop packages
-        run: |
-          julia --project=test -e '
-            using Pkg
-            Pkg.develop([PackageSpec(path=pwd())])
-            Pkg.instantiate()'
-      - uses: julia-actions/julia-buildpkg@v1
+          show-versioninfo: true
+      - uses: julia-actions/cache@v2
+      - uses: julia-actions/julia-runtest@v1
         with:
-          project: test
-      - name: Print package status
-        run: |
-          julia --project=test -e '
-            using InteractiveUtils
-            versioninfo()
-            using Pkg
-            Pkg.status(;mode=Pkg.PKGMODE_MANIFEST)'
-      - name: Run tests
-        run: |
-          cd test && 
-          julia --project --color=yes --depwarn=yes --warn-overwrite=yes --warn-scope=yes --check-bounds=yes --compiled-modules=no  runtests.jl 
+          coverage: false
+          check_bounds: "yes"

.gitignore

@@ -3,6 +3,8 @@
 
 #VS Code
 /.vscode/
+/agent-depot/
+/agent-logs/
 
 *.jl.*.cov
 *.jl.cov

Project.toml

@@ -1,7 +1,7 @@
 name = "SimplexCellLists"
 uuid = "fbe69cd6-f244-4ce4-a90b-8979d86c9ea4"
 authors = ["Nathan Zimmerberg"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
@@ -10,5 +10,9 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 ArgCheck = "2"
+LinearAlgebra = "1"
 StaticArrays = "1"
-julia = "1.8"
+julia = "1.12"
+
+[workspace]
+projects = ["test", "benchmark"]

README.md

@@ -1,162 +1,31 @@
 # SimplexCellLists WIP
 
 [![Build Status](https://github.com/medyan-dev/SimplexCellLists.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/medyan-dev/SimplexCellLists.jl/actions/workflows/CI.yml?query=branch%3Amain)
-![Julia Support](https://img.shields.io/badge/Julia-1.8-blue)
 
-This Julia package accelerates computations on all pairs of 3D points, line segments, and triangles within a cutoff distance.
+This Julia package contains data structures and algorithms for doing computations on pairs of 3D points, line segments, and triangles within a cutoff distance.
 
-This package is largely inspired by [CellListMap.jl](https://github.com/m3g/CellListMap.jl).
-
-However, there is no support for periodic boundary conditions, 2D systems, or types other than Float32.
-
-See [CellListMap.jl](https://github.com/m3g/CellListMap.jl) if you want these features, or higher performance on systems without triangles and line segments.
-
-
-## Setup
-
-- `Simplex{N}` is `SVector{N, SVector{3, Float32}}`
-- `Point` is `Simplex{1}`
-- `Line` is `Simplex{2}`
-- `Triangle` is `Simplex{3}`
-
-There are multiple algorithms that implement the `SimplexCellList` abstract type interface.
-Currently:
-- `Naive`: A simple `SimplexCellList` with no extra options that uses double for loops.
-- `Painter`: A grid is painted with element ids based on a max range. Based on the ideas in [cytosim](https://gitlab.com/f-nedelec/cytosim/-/blob/af739d2ff768628e4737d3a75457676e1a7f4287/src/sim/fiber_grid.h).
-
-Let `T` be a concrete subtype of `SimplexCellList`
-
-### Constructor
-Construct `T`:
-
-```julia
-T(numpointgroups::Integer, numlinegroups::Integer, numtrianglegroups::Integer; kwargs...)::T
-```
-
-Where `numpointgroups` is the number of groups of points, `numlinegroups` is the number of groups of lines,
-and `numtrianglegroups` is the number of groups of triangles.
-
-`kwargs` are options specific for `T`
-
-### `setElements!`
-Reset the elements stored in `s` in batch:
-
-```julia
-setElements!(s::T, points, lines, triangles)::Nothing
-```
-
-Where `points`, `lines` and `triangles` are collections of collections of objects convertible to 
-`Point`, `Line`, and `Triangle` respectively.
-
-For example, each collection in `points` is a group of points that can be mapped over independently from, or together with, other groups.
-
-Added elements will have a group index and element index based on the order of the inputs.
-The first group in each type has group index 1, and the first element in each group has element index 1.
-
-### `addElement!`
-Add a new element to `s`, and return its element index:
-
-```julia
-addElement!(s::T, group_idx::Integer, element::Simplex{N})::Int32
-```
-The new element will be pushed to the end of the specified group.
-
-### `deactivate!`
-Deactivate an existing element in `s`
-
-```julia
-deactivate!(s::T, group_idx::Integer, element_idx::Integer, element_type::Type{Simplex{N}})::Nothing
-```
-Inactive elements are not mapped over. Elements are active by default.
-
-### `activate!`
-Re-activate an existing element in `s`
+It provides cell lists for points ([`PointCellList`](src/pointcelllist.jl)) and line segments ([`LineSegCellList`](src/linesegcelllist.jl)) with fast nearby-neighbor mapping, and squared-distance functions (`dist_sqr`) for all pairs of points, line segments, and triangles.
 
-```julia
-activate!(s::T, group_idx::Integer, element_idx::Integer, element_type::Type{Simplex{N}})::Nothing
-```
-Inactive elements are not mapped over. Elements are active by default.
-
-### `isActive`
-Return if an existing element in `s` is active.
-
-```julia
-isActive(s::T, group_idx::Integer, element_idx::Integer, element_type::Type{Simplex{N}})::Bool
-```
-Inactive elements are not mapped over. Elements are active by default.
-
-## Mapping
-
-The following functions allow mapping a custom function over pairs of simplexes within some cutoff.
-
-### Mapped function `f`
-
-The function f should have the same form as used in CellListMap.jl.
-
-`i` is the element index of simplex `x`, `j` is the element index of simplex `y`. 
-
-`d2` is an approximate `Float32` squared distance between `x` and `y`.
-
-Except here `x` and `y` are `Simplex{N}`, `Simplex{M}`
-
-```julia
-    function f(x,y,i,j,d2,output)
-        # update output
-        return output
-    end
-```
-
-The order in which pairs of elements in range are mapped is implementation dependent.
-
-The elements passed to `f` may be slightly different from the elements added to `T` due to implementation dependent floating point rounding errors.
-
-If a pair distance is very near the cutoff, it is implementation dependent whether the pair gets mapped or not due to floating point rounding errors.
-
-Therefore, if more precision is needed, add some extra distance to the cutoff, store the elements externally in 64 bit precision, and in `f` use `i` and `j` to get the precise elements and again check distances.
-
-### `mapSimplexElements`
-
-Map `f` to all simplexes in a group close to a single simplex.
-
-```julia
-mapSimplexElements(f, output, s::T, x::Simplex{N}, group_idx::Integer, elements_type::Type{Simplex{M}}, cutoff::Float32) where {N, M}
-```
-
-Apply function `f` to all elements in group `group_idx` within the cutoff range of the simplex `x`, and
-return the output of the final `f` call.
+This package is largely inspired by [CellListMap.jl](https://github.com/m3g/CellListMap.jl).
 
-`x` is always `x` and `i` is always 0, in calls to `f`.
+However, there is no support for periodic boundary conditions or 2D systems.
 
-### `mapPairElements`
+See [CellListMap.jl](https://github.com/m3g/CellListMap.jl) if you want these features.
 
-Map `f` to all pairs of nearby simplexes in a single group.
+## Quickstart
 
 ```julia
-mapPairElements(f, output, s::T, group_idx::Integer, elements_type::Type{Simplex{N}}, cutoff::Float32) where {N}
-```
-Apply function `f` to all unordered pairs of elements in group `group_idx` within cutoff range, and return the output of the final `f` call.
-
-`f` is never called more than once per unordered pair. Which element is `x` and `y` in calls to `f` is implementation dependent.
-
+using SimplexCellLists, StaticArrays
 
-### `mapElementsElements`
+# 10×10×10 grid of cells, each 1.0 wide, centered on the origin
+pcl = PointCellList{Int64,Float32}((10, 10, 10), 1.0f0)
+cell_point_add!(pcl, SA[0.1f0, 0.2f0, 0.3f0], 1)
+cell_point_add!(pcl, SA[0.3f0, 0.2f0, 0.3f0], 2)
 
-Map `f` to all pairs of nearby simplexes between two different groups.
-
-```julia
-mapElementsElements(
-        f, 
-        output, 
-        s::T, 
-        x_group_idx::Integer, 
-        x_type::Type{Simplex{N}}, 
-        y_group_idx::Integer, 
-        y_type::Type{Simplex{M}}, 
-        cutoff::Float32,
-    ) where {N, M}
+# Count the points within 0.5 of the origin
+n = map_nearby_points(pcl, SA[0.0f0, 0.0f0, 0.0f0], 0.5f0, 0) do entry, sep, count
+    (count + 1, true)
+end
 ```
-Apply function `f` to each pair of elements from two different groups that are within cutoff range of each other, and return the output of the final `f` call.
-
-The first element has is an `x_type` in group `x_group_idx` and the second element is a `y_type` in group `y_group_idx`.
 
-`f` is never called more than once per pair.
+`LineSegCellList`, `cell_line_seg_add!`, and `map_nearby_line_segs` are the line segment analogs, and `cell_points_clear!`/`cell_line_segs_clear!` reset a list for reuse.

benchmark/Project.toml

@@ -0,0 +1,6 @@
+[deps]
+Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SimplexCellLists = "fbe69cd6-f244-4ce4-a90b-8979d86c9ea4"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

benchmark/distancesbench.jl

@@ -2,11 +2,11 @@ using StaticArrays
 using LinearAlgebra
 using SimplexCellLists
 using Random
-using BenchmarkTools
+using Chairmarks
 
 function dist2linelinevect!(r,a,b)
     @inbounds for i in eachindex(r, a, b)
-        @inline r[i] = SimplexCellLists.dist2LineLine(a[i],b[i])
+        @inline r[i] = dist_sqr(a[i],b[i])
     end
 end
 
@@ -15,13 +15,13 @@ a = rand(SVector{2,SVector{3,Float32}}, N)
 b = rand(SVector{2,SVector{3,Float32}}, N)
 r = zeros(Float32, N)
 println("line line")
-@btime dist2linelinevect!($r,$a,$b)
+display(@b dist2linelinevect!($r,$a,$b))
 println()
 
 
 function dist2pointlinevect!(r,a,b)
     @inbounds for i in eachindex(r, a, b)
-        @inline r[i] = SimplexCellLists.dist2PointLine(a[i],b[i])
+        @inline r[i] = dist_sqr(a[i],b[i])
     end
 end
 
@@ -30,12 +30,12 @@ a = rand(SVector{1,SVector{3,Float32}}, N)
 b = rand(SVector{2,SVector{3,Float32}}, N)
 r = zeros(Float32, N)
 println("point line")
-@btime dist2pointlinevect!($r,$a,$b)
+display(@b dist2pointlinevect!($r,$a,$b))
 println()
 
 function dist2pointpointvect!(r,a,b)
     @inbounds for i in eachindex(r, a, b)
-        @inline r[i] = SimplexCellLists.dist2PointPoint(a[i],b[i])
+        @inline r[i] = dist_sqr(a[i],b[i])
     end
 end
 
@@ -44,12 +44,12 @@ a = rand(SVector{1,SVector{3,Float32}}, N)
 b = rand(SVector{1,SVector{3,Float32}}, N)
 r = zeros(Float32, N)
 println("point point")
-@btime dist2pointpointvect!($r,$a,$b)
+display(@b dist2pointpointvect!($r,$a,$b))
 println()
 
 function dist2pointtrianglevect!(r,a,b)
     @inbounds for i in eachindex(r, a, b)
-        @inline r[i] = SimplexCellLists.dist2PointTriangle(a[i],b[i])
+        @inline r[i] = dist_sqr(a[i],b[i])
     end
 end
 
@@ -58,12 +58,12 @@ a = rand(SVector{1,SVector{3,Float32}}, N)
 b = rand(SVector{3,SVector{3,Float32}}, N)
 r = zeros(Float32, N)
 println("point triangle")
-@btime dist2pointtrianglevect!($r,$a,$b)
+display(@b dist2pointtrianglevect!($r,$a,$b))
 println()
 
 function dist2linetrianglevect!(r,a,b)
     @inbounds for i in eachindex(r, a, b)
-        @inline r[i] = SimplexCellLists.dist2LineTriangle(a[i],b[i])
+        @inline r[i] = dist_sqr(a[i],b[i])
     end
 end
 
@@ -72,12 +72,12 @@ a = rand(SVector{2,SVector{3,Float32}}, N)
 b = rand(SVector{3,SVector{3,Float32}}, N)
 r = zeros(Float32, N)
 println("line triangle")
-@btime dist2linetrianglevect!($r,$a,$b)
+display(@b dist2linetrianglevect!($r,$a,$b))
 println()
 
 function dist2triangletrianglevect!(r,a,b)
     @inbounds for i in eachindex(r, a, b)
-        @inline r[i] = SimplexCellLists.dist2TriangleTriangle(a[i],b[i])
+        @inline r[i] = dist_sqr(a[i],b[i])
     end
 end
 
@@ -86,5 +86,5 @@ a = rand(SVector{3,SVector{3,Float32}}, N)
 b = rand(SVector{3,SVector{3,Float32}}, N)
 r = zeros(Float32, N)
 println("triangle triangle")
-@btime dist2triangletrianglevect!($r,$a,$b)
+display(@b dist2triangletrianglevect!($r,$a,$b))
 println()