-
Notifications
You must be signed in to change notification settings - Fork 0
/
dispersal.R
187 lines (184 loc) · 7 KB
/
dispersal.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
# Copyright (C) 2023, 2024 Stefan Fallert, Lea Li, Juliano Sarmento Cabral
#
# This file is part of metaRange.
#
# metaRange is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 3.
#
# metaRange is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with metaRange. If not, see <http://www.gnu.org/licenses/>.
#' Negative Exponential kernel
#'
#' @param x `<numeric>` distance at which the probability should be calculated.
#' @param mean_dispersal_dist `<numeric>` mean dispersal distance. Needs to be (>0).
#' @details
#' The negative exponential kernel is defined as:
#' \deqn{f(x) = \frac{1}{2 \pi a^2} e^{-\frac{x}{a}}}{fx = 1 / (2 * pi * a^2) * exp(-x / a)}
#' where \eqn{a} is the mean dispersal distance divided by 2.
#' @examples
#' negative_exponential_function(1, 1)
#' @references
#' Nathan, R., Klein, E., Robledo-Arnuncio, J.J. and Revilla, E. (2012)
#' Dispersal kernels: review.
#' in: *Dispersal Ecology and Evolution* pp. 187--210.
#' (eds J. Clobert, M. Baguette, T.G. Benton and J.M. Bullock),
#' Oxford, UK: Oxford Academic, 2013.
#' \doi{10.1093/acprof:oso/9780199608898.003.0015}
#' @return `<numeric>` The probability at distance x.
#' @export
negative_exponential_function <- function(x, mean_dispersal_dist) {
a <- mean_dispersal_dist / 2
(1 / (2 * pi * (a^2))) * exp(-x / a)
}
#' Calculate 2D dispersal kernel.
#'
#' Use a probability function to create a 2D dispersal kernel matrix.
#'
#' @param max_dispersal_dist `<numeric>` maximum dispersal distance in grid cells.
#' The size (rows and columns) of the created dispersal kernel matrix will be
#' `2 * max_dispersal_dist + 1`.
#' @param kfun `<function>` the probability function that is used to calculate
#' the values for each cell of the dispersal kernel. Can be user-defined,
#' in which case it needs to vectorized and accept (at least) the parameter
#' "x" representing the distance from the source as its input and return a
#' vector of the same size as `max_dispersal_dist`.
#' @param normalize `<boolean>` whether to normalize the kernel (`sum(kernel) == 1)`).
#' @param ... additional parameters to be passed to the kernel function.
#' @details This function first creates an matrix of size `2 * max_dispersal_dist + 1`,
#' where each cell contains the distance from the center of the cell to the center
#' of the matrix. After that, the kernel function (`kfun`) is called with this
#' matrix as input and expected to return a vector with the calculated
#' probabilities for each cell. Lastly, the kernel is optionally normalized.
#' @examples
#' # a very simple uniform kernel
#' uniform_kernel <- calculate_dispersal_kernel(
#' max_dispersal_dist = 2,
#' kfun = function(x) {
#' rep(1, length(x))
#' },
#' normalize = FALSE
#' )
#' # same as
#' stopifnot(
#' uniform_kernel == matrix(1, nrow = 5, ncol = 5)
#' )
#'
#' # How does the input matrix look like,
#' # that is used as input for `kfun`?
#' calculate_dispersal_kernel(
#' max_dispersal_dist = 2,
#' kfun = function(x) {
#' return(x)
#' },
#' normalize = FALSE
#' )
#'
#' # now a negative exponential kernel.
#' # note that `mean_dispersal_dist` is a parameter of
#' # the `negative_exponential_function`
#' # and is passed to `kfun` via the `...`.
#' calculate_dispersal_kernel(
#' max_dispersal_dist = 2,
#' kfun = negative_exponential_function,
#' mean_dispersal_dist = 1
#' )
#' @return Dispersal kernel with probabilities.
#' @export
calculate_dispersal_kernel <- function(
max_dispersal_dist,
kfun,
normalize = TRUE,
...) {
max_dispersal_dist <- checkmate::assert_int(max_dispersal_dist, coerce = TRUE, lower = 1L)
checkmate::assert_function(kfun, args = c("x"))
checkmate::assert_flag(normalize)
size <- 2 * max_dispersal_dist + 1
dispersal_kernel <- matrix(0, nrow = size, ncol = size)
midpoint <- max_dispersal_dist + 1
for (i in 1:size) {
for (j in 1:size) {
dispersal_kernel[i, j] <- sqrt(abs(i - midpoint)^2 + abs(j - midpoint)^2)
}
}
dispersal_kernel <- kfun(x = dispersal_kernel, ...)
if (length(dispersal_kernel) != size^2) {
stop("Kernel function should return: ", size^2, " values, but returned: ", length(dispersal_kernel), ".")
}
dim(dispersal_kernel) <- c(size, size)
if (normalize) {
dispersal_kernel <- dispersal_kernel / sum(dispersal_kernel)
}
return(dispersal_kernel)
}
#' Dispersal process
#'
#' Disperse a (abundance) matrix using a dispersal kernel and optional weights.
#' @param dispersal_kernel `<numeric matrix>` dispersal kernel. A 2D matrix of
#' uneven size, containing the weights that deciedes how the individuals from the
#' cell in the center are going to be distributed to the sourrounding cells.
#' @param abundance `<numeric matrix>` abundance matrix.
#' @param weights `<numeric matrix>` optional weights in form of a matrix
#' that has the same dimensions as the abundance and a range between `0` and `1`.
#' Should not contain any `NA`.
#' @details
#' Each cell in the abundance matrix is dispersed using the dispersal kernel.
#' If a matrix of weights is supplied, the individuals will redistribute
#' within the dispersal kernel according to the weights.
#' I.e. individuals will more likely move towards areas with a higher
#' weight, if they are within their dispersal distance.
#' Note:
#' * the abundance is modified in place, to optimize performance.
#' * Any `NA` or `NaN` in abundance or weights will be (in-place) replaced by `0`.
#' @examples
#' n <- 10
#' n2 <- n^2
#' abu <- matrix(1:n2, nrow = n, ncol = n)
#' suitab <- matrix(1, nrow = n, ncol = n)
#' kernel <- calculate_dispersal_kernel(
#' max_dispersal_dist = 4,
#' kfun = negative_exponential_function,
#' mean_dispersal_dist = 1.2
#' )
#' res1 <- dispersal(
#' dispersal_kernel = kernel,
#' abundance = abu
#' )
#' res2 <- dispersal(
#' dispersal_kernel = kernel,
#' abundance = abu,
#' weights = suitab
#' )
#' stopifnot(sum(res1) - sum(res2) < 0.01)
#' # Note that the abundance is modified in place, i.e:
#' stopifnot(sum(abu - res2) < 0.01)
#' @return `<numeric matrix>` Dispersed abundance matrix.
#' @export
dispersal <- function(
dispersal_kernel,
abundance,
weights) {
if (missing(dispersal_kernel)) {
stop("No dispersal kernel supplied.")
}
if (missing(abundance)) {
stop("No abundance matrix supplied.")
}
if (missing(weights)) {
return(dispersal_fixed_unweighted(
dispersal_kernel = dispersal_kernel,
abundance = abundance
))
} else {
return(dispersal_fixed_weighted(
dispersal_kernel = dispersal_kernel,
abundance = abundance,
weights = weights
))
}
}