Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fractal Triad #66

Merged
merged 12 commits into from
May 8, 2024
3 changes: 1 addition & 2 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
13 changes: 2 additions & 11 deletions docs/src/vignettes/entropize.md
Original file line number Diff line number Diff line change
Expand Up @@ -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))
```
27 changes: 9 additions & 18 deletions docs/src/vignettes/overview.md
Original file line number Diff line number Diff line change
Expand Up @@ -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]
```

Expand All @@ -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
Expand All @@ -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()
```
14 changes: 6 additions & 8 deletions docs/src/vignettes/uniqueness.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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()
```
23 changes: 23 additions & 0 deletions ext/SDTExt.jl
Original file line number Diff line number Diff line change
@@ -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
14 changes: 4 additions & 10 deletions src/BiodiversityObservationNetworks.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,9 @@ export AdaptiveSpatial
include("cubesampling.jl")
export CubeSampling

include("fractaltriad.jl")
export FractalTriad

include("uniqueness.jl")
export Uniqueness

Expand All @@ -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
65 changes: 39 additions & 26 deletions src/adaptivespatial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]))
Expand All @@ -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
Expand All @@ -65,7 +57,7 @@ function _generate!(
end
coords[i] = popat!(pool, best_s)
end
return (coords, uncertainty)
return coords
end

function _matérn(d, ρ, ν)
Expand All @@ -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
49 changes: 25 additions & 24 deletions src/balancedacceptance.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

# ====================================================
Expand All @@ -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
Loading
Loading