Merge pull request #66 from PoisotLab/mdc/fractal_triad

Fractal Triad

Project.toml

@@ -24,8 +24,7 @@ SpeciesDistributionToolkit = "72b53823-5c0b-4575-ad0e-8e97227ad13b"
 NeutralLandscapes = "71847384-8354-4223-ac08-659a5128069f"
 
 [extensions]
-SDMToolkitExt = ["SpeciesDistributionToolkit"]
-NLExt = ["NeutralLandscapes"]
+SDTExt = ["SpeciesDistributionToolkit"]
 
 [compat]
 Distributions = "0.25"

docs/src/vignettes/entropize.md

@@ -29,15 +29,6 @@ pixel scale:
 
 ```@example 1
 U = entropize(measurements)
-heatmap(U)
-```
-
-The values closest to the median of the distribution have the highest entropy, and the values closest to its extrema have an entropy of 0. The entropy matrix is guaranteed to have values on the unit interval.
-
-We can use `entropize` as part of a pipeline, and overlay the points optimized based on entropy on the measurement map:
-
-```@example 1
 locations =
-    measurements |> entropize |> seed(BalancedAcceptance(; numpoints = 100)) |> first
-heatmap(U)
-```
+    seed(BalancedAcceptance(; numsites = 100, uncertainty=U))
+```

docs/src/vignettes/overview.md

@@ -31,31 +31,23 @@ less) uncertainty. To start with, we will extract 200 candidate points, *i.e.*
 
 
 ```@example 1
-pack = seed(BalancedAcceptance(; numpoints = 200), U);
+candidates = seed(BalancedAcceptance(; numsites = 200));
 ```
 
-The output of a `BONSampler` (whether at the seeding or refinement step) is
-always a tuple, storing in the first position a vector of `CartesianIndex`
-elements, and in the second position the matrix given as input. We can have a
-look at the first five points: 
+We can have a look at the first five points: 
 
 ```@example 1
-first(pack)[1:5]
+candidates[1:5]
 ```
 
-Although returning the input matrix may seem redundant, it actually allows to
-chain samplers together to build pipelines that take a matrix as input, and
-return a set of places to sample as outputs; an example is given below.
-
 The positions of locations to sample are given as a vector of `CartesianIndex`,
 which are coordinates in the uncertainty matrix. Once we have generated a
 candidate proposal, we can further refine it using a `BONRefiner` -- in this
 case, `AdaptiveSpatial`, which performs adaptive spatial sampling (maximizing
 the distribution of entropy while minimizing spatial auto-correlation).
 
 ```@example 1
-candidates, uncertainty = pack
-locations, _ = refine(candidates, AdaptiveSpatial(; numpoints = 50), uncertainty)
+locations = refine(candidates, AdaptiveSpatial(; numsites = 50, uncertainty=U))
 locations[1:5]
 ```
 
@@ -72,10 +64,8 @@ functions have a curried version that allows chaining them together using pipes
 
 ```@example 1
 locations =
-    U |>
-    seed(BalancedAcceptance(; numpoints = 200)) |>
-    refine(AdaptiveSpatial(; numpoints = 50)) |>
-    first
+    seed(BalancedAcceptance(; numsites = 200)) |>
+    refine(AdaptiveSpatial(; numsites = 50, uncertainty=U))
 ```
 
 This works because `seed` and `refine` have curried versions that can be used
@@ -84,5 +74,6 @@ the original uncertainty matrix:
 
 ```@example 1
 plt = heatmap(U)
-#scatter!(plt, [x[1] for x in locations], [x[2] for x in locations], ms=2.5, mc=:white, label="")
-```
+scatter!(plt, [x[1] for x in locations], [x[2] for x in locations])
+current_figure()
+```

docs/src/vignettes/uniqueness.md

@@ -32,8 +32,6 @@ Now we'll use the `stack` function to combine our four environmental layers into
 layers = BiodiversityObservationNetworks.stack([temp,precip,elevation]);
 ```
 
-this requires NeutralLandscapes v0.1.2
-
 ```@example 1
 uncert = rand(MidpointDisplacement(0.8), size(temp), mask=temp);
 heatmap(uncert) 
@@ -42,15 +40,15 @@ heatmap(uncert)
 Now we'll get a set of candidate points from a BalancedAcceptance seeder that has no bias toward higher uncertainty values.
 
 ```@example 1
-candpts, uncert = uncert |> seed(BalancedAcceptance(numpoints=100, α=0.0)); 
+candpts = seed(BalancedAcceptance(numsites=100)); 
 ```
 
 Now we'll `refine` our `100` candidate points down to the 30 most environmentally unique.
 
 ```@example 1
-finalpts, uncert = refine(candpts, Uniqueness(;numpoints=30, layers=layers), uncert)
-#=
+finalpts = refine(candpts, Uniqueness(;numsites=30, layers=layers))
 heatmap(uncert)
-scatter!([p[2] for p in candpts], [p[1] for p in candpts], fa=0.0, msc=:white, label="Candidate Points")
-scatter!([p[2] for p in finalpts], [p[1] for p in finalpts], c=:dodgerblue, msc=:white, label="Selected Points")=#
-```
+scatter!([p[1] for p in candpts], [p[2] for p in candpts], color=:white)
+scatter!([p[1] for p in finalpts], [p[2] for p in finalpts], color=:dodgerblue, msc=:white)
+current_figure()
+```

ext/SDTExt.jl

@@ -0,0 +1,23 @@
+module SDTExt
+
+using BiodiversityObservationNetworks
+using SpeciesDistributionToolkit
+
+@info "Loading BONs.jl support for SimpleSDMLayers.jl ..."
+
+function stack(
+    layers::Vector{<:SpeciesDistributionToolkit.SimpleSDMLayers.SimpleSDMLayer},
+)
+    # assert all layers are the same size if we add this function to BONs.jl
+    mat = zeros(size(first(layers))..., length(layers))
+    for (l, layer) in enumerate(layers)
+        thismin, thismax = extrema(layer)
+        mat[:, :, l] .= broadcast(
+            x -> isnothing(x) ? NaN : (x - thismin) / (thismax - thismin),
+            layer.grid,
+        )
+    end
+    return mat
+end
+
+end

src/BiodiversityObservationNetworks.jl

@@ -37,6 +37,9 @@ export AdaptiveSpatial
 include("cubesampling.jl")
 export CubeSampling
 
+include("fractaltriad.jl")
+export FractalTriad
+
 include("uniqueness.jl")
 export Uniqueness
 
@@ -49,16 +52,7 @@ export refine, refine!
 include("entropize.jl")
 export entropize, entropize!
 
-include("optimize.jl")
-export optimize
-
 include("utils.jl")
-export stack, squish, _squish, _norm
-
-#=using Requires
-function __init__()
-    @require SpeciesDistributionToolkit="72b53823-5c0b-4575-ad0e-8e97227ad13b" include(joinpath("integrations", "simplesdms.jl"))
-    @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" include(joinpath("integrations", "zygote.jl"))
-end=#
+export stack
 
 end

src/adaptivespatial.jl

@@ -3,40 +3,32 @@
 
 ...
 
-**numpoints**, an Integer (def. 50), specifying the number of points to use.
-
-**α**, an AbstractFloat (def. 1.0), specifying ...
+**numsites**, an Integer (def. 50), specifying the number of points to use.
 """
-Base.@kwdef mutable struct AdaptiveSpatial{T <: Integer} <: BONRefiner
-    numpoints::T = 50
-    function AdaptiveSpatial(numpoints)
-        if numpoints < one(numpoints)
-            throw(
-                ArgumentError(
-                    "You cannot have an AdaptiveSpatial with fewer than one point",
-                ),
-            )
-        end
-        return new{typeof(numpoints)}(numpoints)
+Base.@kwdef mutable struct AdaptiveSpatial{T <: Integer, F <: AbstractFloat} <: BONRefiner
+    numsites::T = 30
+    uncertainty::Array{F, 2} = rand(50, 50)
+    function AdaptiveSpatial(numsites, uncertainty)
+        as = new{typeof(numsites), typeof(uncertainty[begin])}(numsites, uncertainty)
+        check_arguments(as)
+        return as
     end
 end
 
-_generate!(
-    coords::Vector{CartesianIndex},
-    pool::Vector{CartesianIndex},
-    sampler::AdaptiveSpatial,
-    uncertainty) = throw(ArgumentError(
-        "You can only call AdaptiveSpatial with a single layer of type Matrix"))
+maxsites(as::AdaptiveSpatial) = prod(size(as.uncertainty))
+function check_arguments(as::AdaptiveSpatial)
+    check(TooFewSites, as)
+    check(TooManySites, as)
+    return nothing
+end
 
 function _generate!(
     coords::Vector{CartesianIndex},
     pool::Vector{CartesianIndex},
     sampler::AdaptiveSpatial,
-    uncertainty::Array{T,2},
-) where {T <: AbstractFloat}
-
+)
     # Distance matrix (inlined)
-    d = zeros(Float64, Int((sampler.numpoints * (sampler.numpoints - 1)) / 2))
+    d = zeros(Float64, Int((sampler.numsites * (sampler.numsites - 1)) / 2))
 
     # Start with the point with maximum entropy
     imax = last(findmax([uncertainty[i] for i in pool]))
@@ -46,7 +38,7 @@ function _generate!(
     best_score = 0.0
     best_s = 1
 
-    for i in 2:(sampler.numpoints)
+    for i in 2:(sampler.numsites)
         for (ci, cs) in enumerate(pool)
             coords[i] = cs
             # Distance update
@@ -65,7 +57,7 @@ function _generate!(
         end
         coords[i] = popat!(pool, best_s)
     end
-    return (coords, uncertainty)
+    return coords
 end
 
 function _matérn(d, ρ, ν)
@@ -86,3 +78,24 @@ function _D(a1::T, a2::T) where {T <: CartesianIndex{2}}
     x2, y2 = a2.I
     return sqrt((x1 - x2)^2.0 + (y1 - y2)^2.0)
 end
+
+# ====================================================
+#
+#   Tests
+#
+# =====================================================
+
+@testitem "AdaptiveSpatial default constructor works" begin
+    @test typeof(AdaptiveSpatial()) <: AdaptiveSpatial
+end
+
+@testitem "AdaptiveSpatial has right subtypes" begin
+    @test AdaptiveSpatial <: BONRefiner
+    @test AdaptiveSpatial <: BONSampler
+end
+
+@testitem "AdaptiveSpatial requires positive number of sites" begin
+    @test_throws TooFewSites AdaptiveSpatial(numsites = 1)
+    @test_throws TooFewSites AdaptiveSpatial(numsites = 0)
+    @test_throws TooFewSites AdaptiveSpatial(numsites = -1)
+end

src/balancedacceptance.jl

@@ -5,30 +5,34 @@ A `BONSeeder` that uses Balanced-Acceptance Sampling (Van-dem-Bates et al. 2017
 https://doi.org/10.1111/2041-210X.13003)
 """
 Base.@kwdef struct BalancedAcceptance{I <: Integer} <: BONSeeder
-    numpoints::I = 30
-    function BalancedAcceptance(numpoints)
-        bas = new{typeof(numpoints)}(numpoints)
+    numsites::I = 30
+    dims::Tuple{I, I} = (50, 50)
+    function BalancedAcceptance(numsites, dims)
+        bas = new{typeof(numsites)}(numsites, dims)
         check_arguments(bas)
         return bas
     end
 end
 
+maxsites(bas::BalancedAcceptance) = prod(bas.dims)
+
 function check_arguments(bas::BalancedAcceptance)
-    return check(TooFewSites, bas)
+    check(TooFewSites, bas)
+    check(TooManySites, bas)
+    return nothing
 end
 
 function _generate!(
     coords::Vector{CartesianIndex},
-    ::BalancedAcceptance,
-    uncertainty::Matrix{T},
-) where {T <: Real}
+    ba::BalancedAcceptance,
+)
     seed = rand(Int32.(1e0:1e7), 2)
-    x, y = size(uncertainty)
+    x, y = ba.dims
     for (idx, ptct) in enumerate(eachindex(coords))
         i, j = haltonvalue(seed[1] + ptct, 2), haltonvalue(seed[2] + ptct, 3)
         coords[idx] = CartesianIndex(convert.(Int32, [ceil(x * i), ceil(y * j)])...)
     end
-    return (coords, uncertainty)
+    return coords
 end
 
 # ====================================================
@@ -42,39 +46,36 @@ end
 end
 
 @testitem "BalancedAcceptance requires positive number of sites" begin
-    @test_throws TooFewSites BalancedAcceptance(0)
-    @test_throws TooFewSites BalancedAcceptance(1)
+    @test_throws TooFewSites BalancedAcceptance(numsites = 1)
+    @test_throws TooFewSites BalancedAcceptance(numsites = 0)
+    @test_throws TooFewSites BalancedAcceptance(numsites = -1)
 end
 
 @testitem "BalancedAcceptance can't be run with too many sites" begin
     numpts, numcandidates = 26, 25
     @test numpts > numcandidates   # who watches the watchmen?
-    bas = BalancedAcceptance(numpts)
-    dims = Int32.(floor.([sqrt(numcandidates), sqrt(numcandidates)]))
-    uncert = rand(dims...)
-
-    @test_throws TooManySites seed(bas, uncert)
+    dims = Int32.(floor.((sqrt(numcandidates), sqrt(numcandidates))))
+    @test_throws TooManySites BalancedAcceptance(numpts, dims)
 end
 
 @testitem "BalancedAcceptance can generate points" begin
     bas = BalancedAcceptance()
-    sz = (50, 50)
-    coords = seed(bas, rand(sz...)) |> first
+    coords = seed(bas)
 
     @test typeof(coords) <: Vector{CartesianIndex}
-    @test length(coords) == bas.numpoints
+    @test length(coords) == bas.numsites
 end
 
-@testitem "BalancedAcceptance can generate a custom number of points" begin
+@testitem "BalancedAcceptance can generate a custom number of points as positional argument" begin
     numpts = 77
-    bas = BalancedAcceptance(numpts)
     sz = (50, 50)
-    coords = seed(bas, rand(sz...)) |> first
+    bas = BalancedAcceptance(numpts, sz)
+    coords = seed(bas)
     @test numpts == length(coords)
 end
 
 @testitem "BalancedAcceptance can take number of points as keyword argument" begin
     N = 40
-    bas = BalancedAcceptance(; numpoints = N)
-    @test bas.numpoints == N
+    bas = BalancedAcceptance(; numsites = N)
+    @test bas.numsites == N
 end