Merge pull request #57 from EcoJulia/mdc/temporalchange

[new feature] Temporal change with expected statistical properties

.github/workflows/CI.yml

@@ -18,6 +18,8 @@ jobs:
         version:
           - '1.4'
           - '1.5'
+          - '1.6'
+          - '1.7'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,7 +1,7 @@
 name = "NeutralLandscapes"
 uuid = "71847384-8354-4223-ac08-659a5128069f"
 authors = ["Timothée Poisot <timothee.poisot@umontreal.ca>", "Michael Krabbe Borregaard <mkborregaard@sund.ku.dk>", "Michael David Catchen <michael.catchen@mail.mcgill.ca>", "Rafael Schouten <rafaelschouten@gmail.com>", "Virgile Baudrot <virgile.baudrot@posteo.net>"]
-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"

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",

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

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
+
+
+
+

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)],
 )
 
 

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 

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
 

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

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

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

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,21 +71,20 @@ 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
 
     # Mesh
     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)