-
Notifications
You must be signed in to change notification settings - Fork 0
/
get_diversity_timeseries.m
339 lines (286 loc) · 10.8 KB
/
get_diversity_timeseries.m
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
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
%% Make observation in virtual landscape
%
% Returns the diversity timeseries observed within a selected sampling
% grain.
function [div, grain2mean] = get_diversity_timeseries(Stot, c, g, tend, nsamples, change, change_rate_rel, varargin)
%% Preparations & checks
assert(isequal(size(Stot), [1 1]), 'Specify Stot as first input argument.')
% Default: no visualization of sampled grain
visualize = 'off';
if any(strcmp(varargin, 'visualize')) || any(strcmp(varargin, 'plot_grain'))
ind = find(strcmp(varargin, 'visualize'));
visualize = varargin{ind+1};
end
% Default: no plotting of movement
movement = 'off';
if any(strcmp(varargin, 'movement')) || any(strcmp(varargin, 'plot_movement'))
ind = find(strcmp(varargin, 'movement'));
movement = 'on';
end
% Default: no plotting of movement
plot_diversity = 'off';
if any(strcmp(varargin, 'plot_diversity'))
plot_diversity = 'on';
end
% Default: random sampling. Alternatively: site-selection bias.
sampling = 'random';
if any(strcmp(varargin, 'sampling'))
ind = find(strcmp(varargin, 'sampling'));
sampling = varargin{ind+1};
end
% Default: no timelag, all time points used
timelag = 0;
if any(strcmp(varargin, 'timelag'))
ind = find(strcmp(varargin, 'timelag'));
timelag = varargin{ind+1};
end
if timelag == 0
tendL = tend;
else
tendL = tend*timelag;
end
% Default: no truncating. Alternatively: truncate pt first time points of timeseries
truncate = 'off';
pt = 0;
if any(strcmp(varargin, 'truncate'))
ind = find(strcmp(varargin, 'truncate'));
pt = varargin{ind+1};
if pt == 0
truncate = 'off';
else
truncate = 'on';
pt = varargin{ind+1};
end
end
% mean of lognormal SAD distribution
if any(strcmp(varargin, 'mu'))
ind = find(strcmp(varargin, 'mu'));
mu = varargin{ind+1};
end
% variance of lognormal SAD distribution
if any(strcmp(varargin, 'var'))
ind = find(strcmp(varargin, 'var'));
var = varargin{ind+1};
end
% initial spatial distribution of individuals
start_coordinates = 'random';
if any(strcmp(varargin, 'aggregated'))
start_coordinates = 'aggregated';
elseif any(strcmp(varargin, 'heterogeneous'))
start_coordinates = 'heterogeneous';
end
%% Pre-assign empty vectors
% These include only observed time points (minus timelag)
div.S = NaN(nsamples, tend);
div.N = NaN(nsamples, tend);
div.S_PIE = NaN(nsamples, tend);
for s = 1 : nsamples
%% Get a random SAD
% Get species-abundance-distribution (Ni)
% Check whether the variables mu and var exist
if exist('mu', 'var') && exist('var', 'var')
SAD = get_SAD(Stot, 'lognormal', 'mu', mu, 'var', var, 'visualize', 'off');
elseif ~exist('var', 'var') && ~exist('var', 'var')
SAD = get_SAD(Stot, 'lognormal', 'visualize', 'off');
else
error('Did you specify only one of mu and var for SAD?')
end
% for each individual, give species as integer number
species_in = repelem(1:Stot, SAD)';
%% Simulate movement
% get absolute change rate (number of individuals) from relative change
% rate
change_rate = round(change_rate_rel*sum(SAD));
% Simulation
if strcmp(start_coordinates, 'heterogeneous')
% In this case, the SAD changes because it is constructed from 4
% individual SADs
[posXL, posYL, SAD] = simulate_timeseries(SAD, c, g, tendL, change, change_rate, start_coordinates);
species_in = repelem(1:Stot, SAD)';
else
[posXL, posYL] = simulate_timeseries(SAD, c, g, tendL, change, change_rate, start_coordinates);
end
% In the case of "evenness change", the basic movement is default, but
% the change is calculated separately as a changing SAD
if strcmp(change, 'evenness')
species = get_evenness_change(species_in, tend, change_rate);
elseif strcmp(change, 'S')
species = get_richness_change(species_in, tend, change_rate);
else
species = repmat(species_in, 1, tend);
end
%% Get sampling grain
% find lower-left point of all possible gxg rectangles within cxc
% (lower left points range between 0 and c-g)
num_grains = (floor(c/g)^2);
gsteps = (0:((c/g)-1))*g;
[gx_all_mesh, gy_all_mesh] = meshgrid(gsteps, gsteps);
gx_all = gx_all_mesh(:);
gy_all = gy_all_mesh(:);
% get the initial abundance of all grains
comm_t0_N = zeros(num_grains,1);
comm_t0_x = posXL(:,1);
comm_t0_y = posYL(:,1);
for i = 1:length(comm_t0_x)
ind = (gx_all == floor(comm_t0_x(i)/g)*g) & ...
(gy_all == floor(comm_t0_y(i)/g)*g);
comm_t0_N(ind) = comm_t0_N(ind) + 1;
end
% get the initial richness of all grains
visual = get_visual_specifications(species);
rich_t0 = zeros(num_grains,1); % richness per grain
for i = 1:Stot
occs = unique([floor(posXL(species(:,1)==i,1)/g)*g ...
floor(posYL(species(:,1)==i,1)/g)*g], 'rows');
ind = ismember([gx_all gy_all], occs, 'rows');
rich_t0(ind) = rich_t0(ind) + 1;
end
% Choose sampling grain according to sampling strategy
if strcmp(sampling, 'random') % Random grain
rand_idx = randi([1 num_grains], 1, 1);
gx = gx_all(rand_idx);
gy = gy_all(rand_idx);
elseif strcmp(sampling, 'pop-biased') || strcmp(sampling, 'biased') % Population bias
% find the rectangle among these with the highest initial abundance
% of species 1
spec1_t0_x = posXL(species(:,1)==1,1);
spec1_t0_y = posYL(species(:,1)==1,1);
spec1_t0_N = zeros(num_grains,1);
for i = 1:length(spec1_t0_x)
ind = (gx_all == floor(spec1_t0_x(i)/g)*g) & ...
(gy_all == floor(spec1_t0_y(i)/g)*g);
spec1_t0_N(ind) = spec1_t0_N(ind) + 1;
end
[~,idx] = sort(spec1_t0_N, 'descend');
gx = gx_all(idx(1));
gy = gy_all(idx(1));
elseif strcmp(sampling, 'comm-biased') % Community bias
% find the rectangle with the highest initial abundance
% of species (all species, total N in the grain)
[~,idx] = sort(comm_t0_N, 'descend');
gx = gx_all(idx(1));
gy = gy_all(idx(1));
elseif strcmp(sampling, 'rich-biased') % Richness bias
% find the rectangle among these with the highest initial S
[~,idx] = sort(rich_t0, 'descend');
gx = gx_all(idx(1));
gy = gy_all(idx(1));
else
error('Unspecified sampling strategy. Choose "random" or "biased".')
end
%% Save in which quadrant the chosen grain lies
grain2mean.chosen_quadrant(s) = 1*(gx<c/2 && gy>=c/2) ...
+ 2*(gx>=c/2 && gy>=c/2) ...
+ 3*(gx<c/2 && gy<c/2) ...
+ 4*(gx>=c/2 && gy<c/2);
%% Visualize the sampled grain
if strcmp(visualize, 'on')
visual = get_visual_specifications(species);
figure('color', 'white', 'position', [-548,464,342,276])
hold on
set(gca, 'XLim', [0 c], 'YLim', [0 c])
xlabel('Horizontal coordinate')
ylabel('Vertical coordinate')
if strcmp(sampling, 'pop-biased') || strcmp(sampling, 'biased')
splot = scatter(posXL(:,1), posYL(:,1), visual.szs(:,1), visual.colorder(species(:,1),:), 'filled');
splot.MarkerFaceAlpha = 0.5;
scatter(posXL(species(:,1)==1,1), posYL(species(:,1)==1,1), 'r')
title('Sampling grain (red: species 1 at t=0)')
elseif strcmp(sampling, 'random') || strcmp(sampling, 'comm-biased') ...
|| strcmp(sampling, 'rich-biased')
scatter(posXL(:,1), posYL(:,1), visual.szs(:,1), visual.colorder(species(:,1),:), 'filled');
end
plot(gx+[0 g g 0 0], gy+[0 0 g g 0], 'r', 'linewidth', 1.5)
set(gca, 'XTick', 0:g:c, 'YTick', 0:g:c)
grid on
% Cannot change the xticks when grid is on. Plot your own grid!
end
%% Make virtual observation
% If "timelag" is specified, reduce time series respectively (so length
% of time series corresponds to tend in the end)
if timelag > 0
keep = 1:timelag:tendL;
else
keep = 1:tendL;
end
posX = posXL(:,keep);
posY = posYL(:,keep);
% Pre-assign vectors for changing N, Ni,.. in every time step
Stot_max = max(max(species));
Ni = NaN(Stot_max,tend);
N = NaN(1,tend);
pi = NaN(Stot_max,tend);
S = NaN(1,tend);
PIE = NaN(1,tend);
S_PIE = NaN(1,tend);
% Boolean: is the individual within or without grain?
isg = posX>=gx & posX<gx+g & posY>=gy & posY<gy+g;
% Diversity metrics over time
for t = 1 : tend
Ni(1:max(species(:,t)),t) = accumarray(species(:,t), isg(:,t));
Ni(isnan(Ni(:,t)),t) = 0;
N(t) = sum(Ni(:,t));
pi(:,t) = Ni(:,t)/N(t);
S(t) = sum(Ni(:,t)>0);
PIE(t) = 1-sum(pi(:,t).^2);
S_PIE(t) = 1/(1-PIE(t));
end
%%%%%%%%%
% Get true diversity in landscape over time!
%%%%%%%%%
% Boolean: is the individual present at all (or NaN because it will be
% added later due to change)
isp = ~isnan(posX);
% Pre-assign vectors
tNi = NaN(Stot_max,tend);
tN = NaN(1,tend);
tpi = NaN(Stot_max,tend);
tS = NaN(1,tend);
tPIE = NaN(1,tend);
tS_PIE = NaN(1,tend);
% Diversity metrics over time
for t = 1 : tend
tNi(1:max(species(:,t)),t) = accumarray(species(:,t), isp(:,t));
tNi(isnan(tNi(:,t)),t) = 0;
tN(t) = sum(tNi(:,t));
tpi(:,t) = tNi(:,t)/tN(t);
tS(t) = sum(tNi(:,t)>0);
tPIE(t) = 1-sum(tpi(:,t).^2);
tS_PIE(t) = 1/(1-tPIE(t));
end
% plot movement for comparison
if strcmp(movement, 'on')
% Plot the movement with:
plot_movement(SAD, c, g, tend, posX, posY, species, 'wait');
end
div.N(s,:) = N;
div.S(s,:) = S;
div.S_PIE(s,:) = S_PIE;
div.tN(s,:) = tN;
div.tS(s,:) = tS;
div.tS_PIE(s,:) = tS_PIE;
if strcmp(plot_diversity, 'on')
plot_diversity_dynamics(div)
end
%% Get spatial heterogeneity:
% initial (at time = 1) difference of S (or N) in the selected site to
% mean S (or N) in all other sites
% difference to mean abundance
grain2mean.N(s,1) = N(1)-mean(comm_t0_N);
grain2mean.N_st(s,1) = (N(1)-mean(comm_t0_N))/sum(SAD);
% assert(sum(SAD)/num_grains==mean(comm_t0_N), 'not the same, should be the same!')
% difference to mean richness
grain2mean.S(s,1) = S(1)-mean(rich_t0);
grain2mean.S_st(s,1) = (S(1)-mean(rich_t0))/Stot;
end
% Optionally left-truncate the time series (cut first pt percent of
% points)
if strcmp(truncate, 'on')
div.N(:,1:pt) = [];
div.S(:,1:pt) = [];
div.S_PIE(:,1:pt) = [];
div.tN(:,1:pt) = [];
div.tS(:,1:pt) = [];
div.tS_PIE(:,1:pt) = [];
end
end