v0.3

Project.toml

@@ -4,6 +4,7 @@ authors = ["michael catchen <michael.catchen@colorado.edu>"]
 version = "0.2.1"
 
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 HaltonSequences = "13907d55-377f-55d6-a9d6-25ac19e11b95"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
@@ -12,25 +13,27 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SimpleSDMLayers = "2c645270-77db-11e9-22c3-0f302a89c64c"
 SliceMap = "82cb661a-3f19-5665-9e27-df437c7e54c8"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Term = "22787eb5-b846-44ae-b979-8e399b8463ab"
+TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 
 [weakdeps]
-SpeciesDistributionToolkit = "72b53823-5c0b-4575-ad0e-8e97227ad13b"
+NeutralLandscapes = "71847384-8354-4223-ac08-659a5128069f"
 
 [extensions]
-SDMToolkitExt = ["SpeciesDistributionToolkit"]
 
 [compat]
 Distributions = "0.25"
 HaltonSequences = "0.2"
+HiGHS = "1.5"
+JuMP = "1.11"
+ProgressMeter = "1.7.2"
 Requires = "1.3"
 SliceMap = "0.2.7"
 SpecialFunctions = "2.1"
 StatsBase = "0.34"
-ProgressMeter = "1.7.2"
-JuMP = "1.11"
-HiGHS = "1.5"
-julia = "1.9" 
+julia = "1.9"

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,4 @@
+module SDTExt
+
+
+end 

src/BiodiversityObservationNetworks.jl

@@ -1,5 +1,6 @@
 module BiodiversityObservationNetworks
 
+using SimpleSDMLayers
 using Distributions
 using Random
 using HaltonSequences
@@ -10,45 +11,58 @@ using SliceMap
 using JuMP
 using HiGHS
 using LinearAlgebra
+using Term
+using TestItems
 
 include("types.jl")
-export BONSeeder, BONRefiner, BONSampler
+export BONSampler
+export Sites, numsites, pool
+export LayerType, DataLayer, InclusionProbability
+export Layer, Stack
+
+include("sample.jl")
+export sample
+
+include("exceptions.jl")
+export BONException, TooFewSites, TooManySites
+
+include("simplerandom.jl")
+export SimpleRandom
+
+include("spatialstratified.jl")
+export SpatiallyStratified
 
 include("balancedacceptance.jl")
 export BalancedAcceptance
 
-include("adaptivespatial.jl")
-export AdaptiveSpatial
+include("weightedbas.jl")
+export WeightedBalancedAcceptance
+
+include("adaptivehotspot.jl")
+export AdaptiveHotspot
 
 include("cubesampling.jl")
 export CubeSampling
 
+include("fractaltriad.jl")
+export FractalTriad
+
+include("grts.jl")
+export GeneralizedRandomTessellatedStratified
+
 include("uniqueness.jl")
 export Uniqueness
 
+#include("seed.jl")
+#export seed, seed!
 
-include("seed.jl")
-export seed, seed!
-
-include("refine.jl")
-export refine, refine!
+#include("refine.jl")
+#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/adaptivehotspot.jl

@@ -0,0 +1,103 @@
+"""
+    AdaptiveHotspot
+
+- **numsites**: an integer specifying the number of points to use (default is 30)
+- **pool**: the sites that could potentially be picked 
+- **uncertainty**: a `Layer` specifying the current uncertainty at each site
+"""
+Base.@kwdef mutable struct AdaptiveHotspot{T <: Integer, F <: AbstractFloat} <: BONSampler
+    numsites::T = 30
+    uncertainty::Layer = Layer(rand(50, 50))
+    function AdaptiveHotspot(numsites, uncertainty)
+        as = new{typeof(numsites), typeof(uncertainty[begin])}(numsites, uncertainty)
+        check_arguments(as)
+        return as
+    end
+end
+
+AdaptiveHotspot(uncertainty::Matrix{T}; numsites = 30) where T = AdaptiveHotspot(numsites, uncertainty)
+
+maxsites(as::AdaptiveHotspot) = prod(size(as.uncertainty))
+function check_arguments(as::AdaptiveHotspot)
+    check(TooFewSites, as)
+    check(TooManySites, as)
+    return nothing
+end
+
+function _generate!(
+    coords::Vector{CartesianIndex},
+    pool::Vector{CartesianIndex},
+    sampler::AdaptiveHotspot,
+)
+    uncertainty = sampler.uncertainty
+    # Distance matrix (inlined)
+    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]))
+    # Add it to the stack
+    coords[1] = popat!(pool, imax)
+
+    best_score = 0.0
+    best_s = 1
+
+    for i in 2:(sampler.numsites)
+        for (ci, cs) in enumerate(pool)
+            coords[i] = cs
+            # Distance update
+            start_from = Int((i - 1) * (i - 2) / 2) + 1
+            end_at = start_from + Int(i - 2)
+            d_positions = start_from:end_at
+            for ti in 1:(i - 1)
+                d[d_positions[ti]] = _D(cs, coords[ti])
+            end
+            # Get the score
+            score = uncertainty[cs] + sqrt(log(i)) * _h(d[1:end_at], 1.0, 0.5)
+            if score > best_score
+                best_score = score
+                best_s = ci
+            end
+        end
+        coords[i] = popat!(pool, best_s)
+    end
+    return coords
+end
+
+function _matérn(d, ρ, ν)
+    # This is the version from the supp mat
+    # ν = 0.5 to have the exponential version
+    return 1.0 * (2.0^(1.0 - ν)) / gamma(ν) *
+           (sqrt(2ν) * d / ρ)^ν *
+           besselk(ν, sqrt(2ν) * d / ρ)
+end
+
+function _h(d, ρ, ν)
+    K = [_matérn(i, ρ, ν) for i in d]
+    return (0.5 * log(2 * π * ℯ)^length(d)) * sum(K)
+end
+
+function _D(a1::T, a2::T) where {T <: CartesianIndex{2}}
+    x1, y1 = a1.I
+    x2, y2 = a2.I
+    return sqrt((x1 - x2)^2.0 + (y1 - y2)^2.0)
+end
+
+# ====================================================
+#
+#   Tests
+#
+# =====================================================
+
+@testitem "AdaptiveHotspot default constructor works" begin
+    @test typeof(AdaptiveHotspot()) <: AdaptiveHotspot
+end
+
+@testitem "AdaptiveHotspot has right subtypes" begin
+    @test AdaptiveHotspot <: BONSampler
+end
+
+@testitem "AdaptiveHotspot requires positive number of sites" begin
+    @test_throws TooFewSites AdaptiveHotspot(numsites = 1)
+    @test_throws TooFewSites AdaptiveHotspot(numsites = 0)
+    @test_throws TooFewSites AdaptiveHotspot(numsites = -1)
+end