diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index d6eaf94..3851b14 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -18,6 +18,8 @@ jobs: version: - '1.4' - '1.5' + - '1.6' + - '1.7' os: - ubuntu-latest - macOS-latest diff --git a/Project.toml b/Project.toml index e8224ae..97aedec 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NeutralLandscapes" uuid = "71847384-8354-4223-ac08-659a5128069f" authors = ["Timothée Poisot ", "Michael Krabbe Borregaard ", "Michael David Catchen ", "Rafael Schouten ", "Virgile Baudrot "] -version = "0.0.4" +version = "0.1.0" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" @@ -15,7 +15,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] DataStructures = "0.18" Distributions = "0.24, 0.25" -NaNMath = "0.3" +NaNMath = "0.3, 1.0" NearestNeighbors = "0.4" StatsBase = "0.33" julia = "1.4" diff --git a/docs/make.jl b/docs/make.jl index 667738c..8288e91 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,5 +1,8 @@ using Documenter, NeutralLandscapes +# For GR docs bug +ENV["GKSwstype"] = "100" + makedocs( sitename="Neutral Landscapes", authors="Timothée Poisot", diff --git a/docs/src/index.md b/docs/src/index.md index 4baa78e..a16307f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -27,6 +27,22 @@ rand rand! ``` +## Temporal Change +```@docs + NeutralLandscapeUpdater + TemporallyVariableUpdater + SpatiallyAutocorrelatedUpdater + SpatiotemporallyAutocorrelatedUpdater + update + update! + normalize + NeutralLandscapes.rate + NeutralLandscapes.variability + NeutralLandscapes.spatialupdater + NeutralLandscapes._update +``` + + ## Other functions ```@docs diff --git a/src/NeutralLandscapes.jl b/src/NeutralLandscapes.jl index e65d4bd..82c8884 100644 --- a/src/NeutralLandscapes.jl +++ b/src/NeutralLandscapes.jl @@ -1,11 +1,11 @@ module NeutralLandscapes import NaNMath -using StatsBase: sample +using StatsBase: sample, ZScoreTransform, fit, transform using Random: rand! using Statistics: quantile, mean using Distributions: Normal -using NearestNeighbors: KDTree, knn, nn +using NearestNeighbors: KDTree, knn, nn, always_false, knn_point!, SVector using DataStructures: IntDisjointSets, union!, find_root, push! using Base: @kwdef @@ -40,4 +40,20 @@ include("makers/planargradient.jl") include("makers/rectangularcluster.jl") include("makers/wavesurface.jl") +include("updaters/update.jl") +include("updaters/temporal.jl") +include("updaters/spatial.jl") +include("updaters/spatiotemporal.jl") +export update, update! +export NeutralLandscapeUpdater +export TemporallyVariableUpdater +export SpatiallyAutocorrelatedUpdater +export SpatiotemporallyAutocorrelatedUpdater +export rate, variability +export normalize + end # module + + + + diff --git a/src/classify.jl b/src/classify.jl index 70d5d14..6d9d0ea 100644 --- a/src/classify.jl +++ b/src/classify.jl @@ -72,9 +72,9 @@ end blend(clusterarray, array, scaling = 1) = blend(clusterarray, [array], [scaling]) const _neighborhoods = Dict( - :rook => ((1, 0), (0, 1)), - :diagonal => ((1, 0), (0, 1), (1, 1)), - :queen => ((1, 0), (0, 1), (1, 1), (1, -1)), + :rook => [(1, 0), (0, 1)], + :diagonal => [(1, 0), (0, 1), (1, 1)], + :queen => [(1, 0), (0, 1), (1, 1), (1, -1)], ) diff --git a/src/landscape.jl b/src/landscape.jl index f1b4a87..f979e8a 100644 --- a/src/landscape.jl +++ b/src/landscape.jl @@ -43,12 +43,11 @@ replaced by `NaN`. function mask!(array::AbstractArray{<:Float64}, maskarray::AbstractArray{<:Bool}) (size(array) == size(maskarray)) || throw(DimensionMismatch("The dimensions of array, $(size(array)), and maskarray, $(size(maskarray)), must match. ")) array[.!maskarray] .= NaN - array + return array end - # Changes the matrix `mat` so that it is between `0` and `1`. function _rescale!(mat) - mat .-= NaNMath.minimum(mat) - mat ./= NaNMath.maximum(mat) + mn, mx = NaNMath.extrema(mat) + mat .= (mat .- mn) ./ (mx - mn) end diff --git a/src/makers/diamondsquare.jl b/src/makers/diamondsquare.jl index a35506e..a9c7345 100644 --- a/src/makers/diamondsquare.jl +++ b/src/makers/diamondsquare.jl @@ -115,13 +115,13 @@ end # Returns the coordinates for the corners of the subsquare (x,y) given a side-length `sideLength`. function _subsquareCornerCoordinates(x::Int, y::Int, sideLength::Int) - corners = [1 .+ sideLength.*i for i in [(x,y), (x+1, y), (x, y+1), (x+1, y+1)]] + return (1 .+ sideLength.*i for i in ((x,y), (x+1, y), (x, y+1), (x+1, y+1))) end # Runs the diamond step of the `DiamondSquare` algorithm on the square defined by # `corners` on the matrix `mat`. The center of the square is interpolated from the # four corners, and is displaced. The displacement is drawn according to `alg.H` and round using `displace` -function _diamond!(mat, alg, round::Int, corners::AbstractVector{Tuple{Int, Int}}) +function _diamond!(mat, alg, round::Int, corners) centerPt = _centerCoordinate(corners) mat[centerPt...] = _interpolate(mat, corners) + _displace(alg.H, round) end @@ -130,32 +130,32 @@ end # by `corners` on the matrix `mat`. The midpoint of each edge of this square is interpolated # by computing the mean value of the two corners on the edge and the center of the square, and the # displacing it. The displacement is drawn according to `alg.H` and round using `displace` -function _square!(mat, alg::DiamondSquare, round::Int, corners::AbstractVector{Tuple{Int, Int}}) +function _square!(mat, alg::DiamondSquare, round::Int, corners) bottomLeft,bottomRight,topLeft,topRight = corners leftEdge, bottomEdge, topEdge, rightEdge = _edgeMidpointCoordinates(corners) centerPoint = _centerCoordinate(corners) - mat[leftEdge...] = _interpolate(mat, [topLeft,bottomLeft,centerPoint]) + _displace(alg.H, round) - mat[bottomEdge...] = _interpolate(mat, [bottomLeft,bottomRight,centerPoint]) + _displace(alg.H, round) - mat[topEdge...] = _interpolate(mat, [topLeft,topRight,centerPoint]) + _displace(alg.H, round) - mat[rightEdge...] = _interpolate(mat, [topRight,bottomRight,centerPoint]) + _displace(alg.H, round) + mat[leftEdge...] = _interpolate(mat, (topLeft,bottomLeft,centerPoint)) + _displace(alg.H, round) + mat[bottomEdge...] = _interpolate(mat, (bottomLeft,bottomRight,centerPoint)) + _displace(alg.H, round) + mat[topEdge...] = _interpolate(mat, (topLeft,topRight,centerPoint)) + _displace(alg.H, round) + mat[rightEdge...] = _interpolate(mat, (topRight,bottomRight,centerPoint)) + _displace(alg.H, round) end # Runs the square step of the `MidpointDisplacement` algorithm on the square defined # by `corners` on the matrix `mat`. The midpoint of each edge of this square is interpolated # by computing the mean value of the two corners on the edge and the center of the square, and the # displacing it. The displacement is drawn according to `alg.H` and round using `displace` -function _square!(mat, alg::MidpointDisplacement, round::Int, corners::AbstractVector{Tuple{Int, Int}}) +function _square!(mat, alg::MidpointDisplacement, round::Int, corners) bottomLeft,bottomRight,topLeft,topRight = corners leftEdge, bottomEdge, topEdge, rightEdge = _edgeMidpointCoordinates(corners) - mat[leftEdge...] = _interpolate(mat, [topLeft,bottomLeft]) + _displace(alg.H, round) - mat[bottomEdge...] = _interpolate(mat, [bottomLeft,bottomRight]) + _displace(alg.H, round) - mat[topEdge...] = _interpolate(mat, [topLeft,topRight]) + _displace(alg.H, round) - mat[rightEdge...] = _interpolate(mat, [topRight,bottomRight]) + _displace(alg.H, round) + mat[leftEdge...] = _interpolate(mat, (topLeft,bottomLeft)) + _displace(alg.H, round) + mat[bottomEdge...] = _interpolate(mat, (bottomLeft,bottomRight)) + _displace(alg.H, round) + mat[topEdge...] = _interpolate(mat, (topLeft,topRight)) + _displace(alg.H, round) + mat[rightEdge...] = _interpolate(mat, (topRight,bottomRight)) + _displace(alg.H, round) end # Computes the mean of a set of points, represented as a list of indicies to a matrix `mat`. -function _interpolate(mat, points::AbstractVector{Tuple{Int,Int}}) +function _interpolate(mat, points) return mean(mat[pt...] for pt in points) end @@ -169,12 +169,12 @@ end # move from 1.0 to 0 as `round` increases. function _displace(H::Float64, round::Int) σ = 0.5^(round*H) - return(rand(Normal(0, σ))) + return rand(Normal(0, σ)) end # Returns the center coordinate for a square defined by `corners` for the # `DiamondSquare` algorithm. -function _centerCoordinate(corners::AbstractVector{Tuple{Int,Int}}) +function _centerCoordinate(corners) bottomLeft,bottomRight,topLeft,topRight = corners centerX::Int = (_xcoord(bottomLeft)+_xcoord(bottomRight)) ÷ 2 centerY::Int = (_ycoord(topRight)+_ycoord(bottomRight)) ÷ 2 @@ -187,7 +187,7 @@ _xcoord(pt::Tuple{Int,Int}) = pt[1] _ycoord(pt::Tuple{Int,Int}) = pt[2] # Returns an array of midpoints for a square defined by `corners` for the `DiamondSquare` algorithm. -function _edgeMidpointCoordinates(corners::AbstractVector{Tuple{Int,Int}}) +function _edgeMidpointCoordinates(corners) # bottom left, bottom right, top left, top right bottomLeft,bottomRight,topLeft,topRight = corners @@ -196,7 +196,7 @@ function _edgeMidpointCoordinates(corners::AbstractVector{Tuple{Int,Int}}) topEdgeMidpoint::Tuple{Int,Int} = ( (_xcoord(topLeft)+_xcoord(topRight))÷2, _ycoord(topLeft)) rightEdgeMidpoint::Tuple{Int,Int} = ( _xcoord(bottomRight), (_ycoord(bottomRight)+_ycoord(topRight))÷2) - edgeMidpoints = [leftEdgeMidpoint, bottomEdgeMidpoint, topEdgeMidpoint, rightEdgeMidpoint] + edgeMidpoints = (leftEdgeMidpoint, bottomEdgeMidpoint, topEdgeMidpoint, rightEdgeMidpoint) return edgeMidpoints end diff --git a/src/makers/distancegradient.jl b/src/makers/distancegradient.jl index d2e2b1d..c56b217 100644 --- a/src/makers/distancegradient.jl +++ b/src/makers/distancegradient.jl @@ -14,11 +14,24 @@ end function _landscape!(mat, alg::DistanceGradient) @assert maximum(alg.sources) <= length(mat) @assert minimum(alg.sources) > 0 - coordinates = Matrix{Float64}(undef, (2, prod(size(mat)))) - for (i, p) in enumerate(Iterators.product(axes(mat)...)) - coordinates[1:2, i] .= p + coordinates = CartesianIndices(mat) + source_coordinates = map(alg.sources) do c + SVector(Tuple(coordinates[c])) + end + idx = Vector{Int}(undef, 1) + dist = Vector{Float64}(undef, 1) + tree = KDTree(source_coordinates) + _write_knn!(mat, dist, idx, tree, coordinates) + return mat +end + +# Function barrier, somehow we lose type stability with this above +function _write_knn!(mat, dist, idx, tree, coordinates) + sortres = false + for i in eachindex(mat) + point = SVector(Tuple(coordinates[i])) + knn_point!(tree, point, sortres, dist, idx, always_false) + mat[i] = dist[1] end - tree = KDTree(coordinates[:,alg.sources]) - mat[:] .= nn(tree, coordinates)[2] return mat end diff --git a/src/makers/nncluster.jl b/src/makers/nncluster.jl index a52af29..166cec0 100644 --- a/src/makers/nncluster.jl +++ b/src/makers/nncluster.jl @@ -23,11 +23,21 @@ function _landscape!(mat, alg::NearestNeighborCluster) classify!(mat, [alg.p, 1 - alg.p]) replace!(mat, 2.0 => NaN) clusters, nClusters = label(mat, alg.n) - coordinates = _coordinatematrix(clusters) + coordinates = CartesianIndices(clusters) sources = findall(!isnan, vec(clusters)) - tree = KDTree(coordinates[:,sources]) - clusters[:] .= clusters[sources[nn(tree, coordinates)[1]]] + cluster_coordinates = map(sources) do c + SVector(Tuple(coordinates[c])) + end + idx = Vector{Int}(undef, 1) + dist = Vector{Float64}(undef, 1) + tree = KDTree(cluster_coordinates) randvals = rand(nClusters) - mat .= randvals[Int.(clusters)] - mat + sortres = false + for i in eachindex(mat) + point = SVector(Tuple(coordinates[i])) + knn_point!(tree, point, sortres, dist, idx, always_false) + cluster = clusters[sources[idx[1]]] + mat[i] = randvals[Int(cluster)] + end + return mat end diff --git a/src/makers/nnelement.jl b/src/makers/nnelement.jl index a5ba8cf..9204481 100644 --- a/src/makers/nnelement.jl +++ b/src/makers/nnelement.jl @@ -18,26 +18,29 @@ and `k` neighbors. The default is to use three cluster and a single neighbor. end end -function _coordinatematrix(mat) - coordinates = Matrix{Float64}(undef, (2, prod(size(mat)))) - for (i, p) in enumerate(Iterators.product(axes(mat)...)) - coordinates[1:2, i] .= p - end - coordinates -end - function _landscape!(mat, alg::NearestNeighborElement) - fill!(mat, 0.0) clusters = sample(eachindex(mat), alg.n; replace=false) - for (i,n) in enumerate(clusters) - mat[n] = i - end - coordinates = _coordinatematrix(mat) - tree = KDTree(coordinates[:,clusters]) - if alg.k == 1 - mat[:] .= nn(tree, coordinates)[1] + # Preallocate for NearestNeighbors + idx = Vector{Int}(undef, alg.k) + dist = Vector{Float64}(undef, alg.k) + coordinates = CartesianIndices(mat) + cluster_coordinates = map(clusters) do c + SVector(Tuple(coordinates[c])) + end + tree = KDTree(cluster_coordinates) + sortres = false + if alg.k == 1 + for i in eachindex(mat) + point = SVector(Tuple(coordinates[i])) + knn_point!(tree, point, sortres, dist, idx, always_false) + mat[i] = idx[1] + end else - mat[:] .= map(mean, knn(tree, coordinates, alg.k)[1]) + for i in eachindex(mat) + point = SVector(Tuple(coordinates[i])) + knn_point!(tree, point, sortres, dist, idx, always_false) + mat[i] = mean(idx) + end end return mat end diff --git a/src/makers/perlinnoise.jl b/src/makers/perlinnoise.jl index 3265682..9ce523a 100644 --- a/src/makers/perlinnoise.jl +++ b/src/makers/perlinnoise.jl @@ -48,10 +48,11 @@ function _landscape!(A, alg::PerlinNoise) # Generate the Perlin noise noise = zeros(eltype(A), (dim, dim)) - meshbuf = Array{eltype(A),3}(undef, dim, dim, 2) + meshbuf1 = Array{eltype(A),2}(undef, dim, dim) + meshbuf2 = Array{eltype(A),2}(undef, dim, dim) nbufs = ntuple(_->Array{eltype(A),2}(undef, dim, dim), 4) for octave in 0:alg.octaves-1 - octave_noise!(noise, meshbuf, nbufs, alg, octave, dim, dim) + octave_noise!(noise, meshbuf1, meshbuf2, nbufs, alg, octave, dim, dim) end # Randomly extract the desired array size @@ -70,7 +71,7 @@ function _landscape!(A, alg::PerlinNoise) end function octave_noise!( - noise, mesh, (n11, n21, n12, n22), alg::PerlinNoise, octave, nrow, ncol + noise, m1, m2, (n11, n21, n12, n22), alg::PerlinNoise, octave, nrow, ncol ) f(t) = @fastmath 6 * t ^ 5 - 15 * t ^ 4 + 10 * t ^ 3 # Wut @@ -78,13 +79,12 @@ function octave_noise!( rp, cp = alg.periods .* alg.lacunarity^(octave) delta = (rp / nrow, cp / ncol) ranges = range(0, rp-delta[1], length=nrow), range(0, cp-delta[2], length=ncol) - _mesh!(mesh, ranges...) - @fastmath mesh .%= 1 + _rem_meshes!(m1, m2, ranges...) # Gradients # This allocates, but the gradients size changes with octave so it needs to # be a very smart in-place `repeat!` or some kind of generator, and the improvement - # may not be than large (~20%). + # may not be that large (~20%). angles = 2pi .* rand(rp + 1, cp + 1) @fastmath gradients = cat(cos.(angles), sin.(angles); dims=3) d = (nrow ÷ rp, ncol ÷ cp) @@ -100,31 +100,28 @@ function octave_noise!( g222 = @view grad[end-nrow+1:end, end-ncol+1:end, 2] # Ramps - m1 = @view mesh[:, :, 1] - m2 = @view mesh[:, :, 2] n11 .= ((m1 .+ m2 ) .* g111 .+ (m1 .+ m2 ) .* g112) n21 .= ((m1 .-1 .+ m2 ) .* g211 .+ (m1 .-1 .+ m2 ) .* g212) n12 .= ((m1 .+ m2 .- 1) .* g121 .+ (m1 .+ m2 .- 1) .* g122) n22 .= ((m1 .-1 .+ m2 .- 1) .* g221 .+ (m1 .-1 .+ m2 .- 1) .* g222) # Interpolation - mesh .= f.(mesh) - t1 = @view mesh[:, :, 1] - t2 = @view mesh[:, :, 2] + m1 .= f.(m1) + m2 .= f.(m2) noise .+= sqrt(2) .* (alg.persistance ^ octave) .* - ((1 .- t2) .* (n11 .* (1 .- t1) .+ t1 .* n21) .+ - t2 .* (n12 .* (1 .- t1) .+ t1 .* n22)) + ((1 .- m2) .* (n11 .* (1 .- m1) .+ m1 .* n21) .+ + m2 .* (n12 .* (1 .- m1) .+ m1 .* n22)) return noise end -function _mesh!(A, x, y) +function _rem_meshes!(m1, m2, x, y) for (i, ival) in enumerate(x), j in 1:length(y) - A[i, j, 1] = ival + @fastmath m1[i, j] = ival % 1 end for i in 1:length(x), (j, jval) in enumerate(y) - A[i, j, 2] = jval + @fastmath m2[i, j] = jval % 1 end - return A + return end function _view_from_square(source, nrow, ncol) diff --git a/src/makers/planargradient.jl b/src/makers/planargradient.jl index 601b06b..0f1dab2 100644 --- a/src/makers/planargradient.jl +++ b/src/makers/planargradient.jl @@ -19,6 +19,6 @@ function _landscape!(mat, alg::PlanarGradient) where {IT <: Integer} eastness = sin(deg2rad(alg.direction)) southness = -1cos(deg2rad(alg.direction)) rows, cols = axes(mat) - mat .= southness .* rows .+ (eastness .* cols)' + mat .= collect(rows) .* southness .+ cols' .* eastness return mat end diff --git a/src/makers/rectangularcluster.jl b/src/makers/rectangularcluster.jl index f2b6ff6..b97eb18 100644 --- a/src/makers/rectangularcluster.jl +++ b/src/makers/rectangularcluster.jl @@ -19,7 +19,7 @@ end function _landscape!(mat, alg::RectangularCluster) mat .= -1.0 - while minimum(mat) == -1.0 + while any(i -> i === -1.0, mat) width, height = rand(alg.minimum:alg.maximum, 2) row = rand(1:(size(mat,1)-(width-1))) col = rand(1:(size(mat,2)-(height-1))) diff --git a/src/updaters/spatial.jl b/src/updaters/spatial.jl new file mode 100644 index 0000000..523bb74 --- /dev/null +++ b/src/updaters/spatial.jl @@ -0,0 +1,32 @@ + +""" + SpatiallyAutocorrelatedUpdater{SU,R,V} + +A `NeutralLandscapeUpdater` that has a prescribed level of +spatial variation (`variability`) and rate of change (`rate`), +and where the spatial distribution of this change is proportional +to a neutral landscape generated with `spatialupdater` at every time +step. + +TODO: make it possible to fix a given spatial updater at each timestep. +""" +@kwdef struct SpatiallyAutocorrelatedUpdater{SU,R,V} <: NeutralLandscapeUpdater + spatialupdater::SU = DiamondSquare(0.5) + rate::R = 0.1 + variability::V = 0.1 +end + +""" + _update(sau::SpatiallyAutocorrelatedUpdater, mat; transform=ZScoreTransform) + +Updates `mat` using spatially autocorrelated change, using the direction, rate, +and spatial updater parameters from `sau`. + +TODO: doesn't necessarily have to be a ZScoreTransform, could be arbitrary +argument +""" +function _update(sau::SpatiallyAutocorrelatedUpdater, mat) + change = rand(spatialupdater(sau), size(mat)) + delta = rate(sau) .+ variability(sau) .* transform(fit(ZScoreTransform, change), change) + mat .+ delta +end diff --git a/src/updaters/spatiotemporal.jl b/src/updaters/spatiotemporal.jl new file mode 100644 index 0000000..3f115d9 --- /dev/null +++ b/src/updaters/spatiotemporal.jl @@ -0,0 +1,35 @@ + +""" + SpatiotemporallyAutocorrelatedUpdater{SU,R,V} + +A `NeutralLandscapeUpdater` that has a prescribed level of +spatial and temporal variation (`variability`) and rate of change (`rate`), +and where the spatial distribution of this change is proportional +to a neutral landscape generated with `spatialupdater` at every time +step. + +TODO: perhaps spatial and temporal should each have their own variability param +""" +@kwdef struct SpatiotemporallyAutocorrelatedUpdater{SU,R,V} <: NeutralLandscapeUpdater + spatialupdater::SU = DiamondSquare(0.1) + rate::R = 0.1 + variability::V = 0.1 +end + + +""" + _update(stau::SpatiotemporallyAutocorrelatedUpdater, mat) + +Updates `mat` using temporally autocorrelated change, using the direction, rate, +and spatialupdater parameters from `stau`. + +TODO: doesn't necessarily have to be a Normal distribution or ZScoreTransform, +could be arbitrary argument +""" +function _update(stau::SpatiotemporallyAutocorrelatedUpdater, mat) + change = rand(spatialupdater(stau), size(mat)) + temporalshift = rand(Normal(0, variability(stau)), size(mat)) + z = transform(fit(ZScoreTransform, change), change) + delta = rate(stau) .+ variability(stau) * z .+ temporalshift + mat .+ delta +end \ No newline at end of file diff --git a/src/updaters/temporal.jl b/src/updaters/temporal.jl new file mode 100644 index 0000000..6f4f8f8 --- /dev/null +++ b/src/updaters/temporal.jl @@ -0,0 +1,28 @@ +""" + TemporallyVariableUpdater{D,S} <: NeutralLandscapeUpdater + +A `NeutralLandscapeUpdater` that has a prescribed level of temporal variation +(`variability`) and rate of change (`rate`), but no spatial correlation in where +change is distributed. +""" +@kwdef struct TemporallyVariableUpdater{D,R,V} <: NeutralLandscapeUpdater + spatialupdater::D = missing + rate::R = 0.1 + variability::V = 0.1 +end + +""" + _update(tvu::TemporallyVariableUpdater, mat) + +Updates `mat` using temporally autocorrelated change, using the direction and +rate parameters from `tvu`. + +TODO: this doesn't have to be a Normal distribution, could be arbitrary +distribution that is continuous and can have mean 0 (or that can be transformed +to have mean 0) +""" +function _update(tvu::TemporallyVariableUpdater, mat) + U = rand(Normal(0, variability(tvu)), size(mat)) + Δ = rate(tvu) .+ U + mat .+ Δ +end diff --git a/src/updaters/update.jl b/src/updaters/update.jl new file mode 100644 index 0000000..0434d92 --- /dev/null +++ b/src/updaters/update.jl @@ -0,0 +1,90 @@ +""" + NeutralLandscapeUpdater + +NeutralLandscapeUpdater is an abstract type for methods for updating a landscape +matrix +""" +abstract type NeutralLandscapeUpdater end + +""" + spatialupdater(up::NeutralLandscapeUpdater) + +All `NeutralLandscapeUpdater`s have a field `rate` which defines the expected +(or mean) change across all cells per timestep. +""" +rate(up::NeutralLandscapeUpdater) = up.rate + +""" + spatialupdater(up::NeutralLandscapeUpdater) + +All `NeutralLandscapeUpdater`'s have a `spatialupdater` field which is either a +`NeutralLandscapeMaker`, or `Missing` (in the case of temporally correlated +updaters). +""" +spatialupdater(up::NeutralLandscapeUpdater) = up.spatialupdater + +""" + variability(up::NeutralLandscapeUpdater) + +Returns the `variability` of a given `NeutralLandscapeUpdater`. The variability +of an updater is how much temporal variation there will be in a generated +time-series of landscapes. +""" +variability(up::NeutralLandscapeUpdater) = up.variability + +""" + normalize(mats::Vector{M}) + +Normalizes a vector of neutral landscapes `mats` such that all values between 0 +and 1. Note that this does not preserve the `rate` parameter for a given +`NeutralLandscapeUpdater`, and instead rescales it proportional to the +difference between the total maximum and total minimum across all `mats`. +""" +function normalize(mats::Vector{M}) where {M<:AbstractMatrix} + mins, maxs = findmin.(mats), findmax.(mats) + totalmin, totalmax = findmin([x[1] for x in mins])[1], findmax([x[1] for x in maxs])[1] + returnmats = copy(mats) + for (i,mat) in enumerate(mats) + returnmats[i] = (mat .- totalmin) ./ (totalmax - totalmin) + end + return returnmats +end + + +""" + update(updater::T, mat) + +Returns one-timestep applied to `mat` based on the `NeutralLandscapeUpdater` +provided (`updater`). +""" +function update(updater::T, mat) where {T<:NeutralLandscapeUpdater} + _update(updater, mat) +end + + +""" + update(updater::T, mat, n::I) + +Returns a sequence of length `n` where the original neutral landscape `mat` is +updated by the `NeutralLandscapeUpdater` `update` for `n` timesteps. +""" +function update(updater::T, mat, n::I) where {T<:NeutralLandscapeUpdater, I<:Integer} + sequence = [zeros(size(mat)) for _ in 1:n] + sequence[begin] .= mat + for i in 2:n + sequence[i] = _update(updater, sequence[i-1]) + end + sequence +end + + +""" + update!(updater::T, mat) + +Updates a landscape `mat` in-place by directly mutating `mat`. +""" +function update!(updater::T, mat) where {T<:NeutralLandscapeUpdater} + mat .= _update(updater, mat) +end + + diff --git a/test/rand.jl b/test/rand.jl index 2b08295..6f8003e 100644 --- a/test/rand.jl +++ b/test/rand.jl @@ -1,5 +1,9 @@ using NeutralLandscapes, Test +@testset "_rescale" begin + @test NeutralLandscapes._rescale!([3.0 -2.0; 0.0 1.0]) == [1.0 0.0; 0.4 0.6] +end + algorithms = ( DiamondSquare(), DiamondSquare(; H=0.2), diff --git a/test/runtests.jl b/test/runtests.jl index 98d5327..dcabf80 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ using Test, SafeTestsets +@time @safetestset "updaters" begin include("updaters.jl") end @time @safetestset "planar gradient" begin include("planar.jl") end @time @safetestset "rand" begin include("rand.jl") end diff --git a/test/updaters.jl b/test/updaters.jl new file mode 100644 index 0000000..7f51733 --- /dev/null +++ b/test/updaters.jl @@ -0,0 +1,41 @@ +using NeutralLandscapes +using Test + +# test it's running +@test TemporallyVariableUpdater() != π +@test SpatiallyAutocorrelatedUpdater() != π +@test SpatiotemporallyAutocorrelatedUpdater() != π + +function testupdaters(model) + updater = model() + + # Test defaults + @test rate(updater) == 0.1 + @test variability(updater) == 0.1 + + # Test kwargs + updater = model(rate = 1.0, variability=0.05, spatialupdater=MidpointDisplacement(0.5)) + @test rate(updater) == 1.0 + @test variability(updater) == 0.05 + @test typeof(updater.spatialupdater) <: NeutralLandscapeMaker + @test updater.spatialupdater == MidpointDisplacement(0.5) + + # Test updating + env = rand(MidpointDisplacement(0.5), 50, 50) + newenv = update(updater, env) + @test env != newenv + + oldenv = deepcopy(env) + update!(updater, env) + @test env != oldenv +end + + + +models = [ + TemporallyVariableUpdater, + SpatiallyAutocorrelatedUpdater, + SpatiotemporallyAutocorrelatedUpdater +] + +testupdaters.(models)