-
Notifications
You must be signed in to change notification settings - Fork 0
/
get_richness_change.m
49 lines (41 loc) · 1.76 KB
/
get_richness_change.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
%% Change species identities to change richness without changing N
% Imposes the biodiversity change regime "richness change" for the
% specified list of species, specified sampling duration and specified rate
% of biodiversity change.
function [species] = get_richness_change(species_in, tend, change_rate)
% total number of individuals (is constant)
I = size(species_in,1);
% Species-identity of individuals over time matrix
species = NaN(I, tend);
species(:,1) = species_in;
% check whether change rate can be realized
% % % if change_rate>0
Ni = hist(species_in, unique(species_in));
available_slots = Ni-1;
rate_ok = sum(available_slots)>=abs(change_rate)*tend;
assert(rate_ok, 'Richness change cannot be realized. Increase number of individuals, decrease change_rate, or decrease tend.')
% "change_rate" individuals of the most common species are turned into new species
for t = 2:tend
Ni = hist(species(:,t-1), unique(species(:,t-1)));
available_slots = Ni-1;
[~, idx] = sort(available_slots, 'ascend');
newspecs = max(species(:,t-1)) + [1:abs(change_rate)];
losers = fliplr(idx);
species(:,t) = species(:, t-1);
placed = 0;
left = abs(change_rate);
i = 1;
while left>0
nplacable = min(left, available_slots(losers(i)));
species(find(species(:,t)==losers(i), nplacable , 'first'),t) = newspecs(1:nplacable);
placed = placed + nplacable;
left = abs(change_rate)-placed;
newspecs(1:nplacable) = [];
i = i+1;
end
end
% for negative change rates: just flip increase time series
if change_rate < 0
species = fliplr(species);
end
assert(all(sum(species>0)==sum(species(:,1)>0)), 'S changes')