-
Notifications
You must be signed in to change notification settings - Fork 0
/
ga.R
252 lines (244 loc) · 10.4 KB
/
ga.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
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
#' Run genetic algorithms to find stability
#'
#' Outer function for investigating the potential for a complex system to be
#' stabilised by component response rates. For a stretch of values of component
#' number (i.e., system size, S), the function runs some number of iterations,
#' and in each iteration, if a system is not initially found to be stable with
#' a set of component response times, then a genetic algorithm is called in
#' attempt to find a stable solution.
#'
#'@return A table showing frequency of stable complex systems found
#'@param max_sp The maximum system component size (S) to simulate
#'@param iters The number of iterations (complex systems) tried for each S
#'@param int_type The type of interaction, including random (0), competitive
#'(1), mutualist (2), or predator-prey (3)
#'@param rmx The standard devation of interaction (off-diagonal) elements in a
#'complex system
#'@param C The connectedness of a complex system (i.e., the probability that an
#'off-diagonal element does not equal zero)
#'@examples
#'sim <- Evo_rand_gen_var(max_sp = 2, iters = 1);
#'@export
Evo_rand_gen_var <- function(max_sp, iters, int_type = 0, rmx = 0.4, C = 1){
tot_res <- NULL;
fea_res <- NULL;
for(i in 2:max_sp){
nn <- i;
A1_stt <- 0;
A2_stt <- 0;
A1_fet <- 0;
A2_fet <- 0;
iter <- iters;
tot_res[[i-1]] <- matrix(data = 0, nrow = iter, ncol = 3);
fea_res[[i-1]] <- matrix(data = 0, nrow = iter, ncol = 2);
while(iter > 0){
r_vec <- rnorm(n = i, mean = 0, sd = rmx);
A0_dat <- rnorm(n = i * i, mean = 0, sd = 0.4);
A0 <- matrix(data = A0_dat, nrow = i, ncol = i);
A0 <- species_interactions(mat = A0, type = int_type);
C_dat <- rbinom(n = i * i, size = 1, prob = C);
C_mat <- matrix(data = C_dat, nrow = i, ncol = i);
A0 <- A0 * C_mat;
diag(A0) <- -1;
gam1 <- runif(n = i, min = 0, max = 2);
A1 <- A0 * gam1;
A0_stb <- max(Re(eigen(A0)$values)) < 0;
A1_stb <- rand_mat_ga(A1);
A0_fea <- min(-1*solve(A0) %*% r_vec) > 0;
A1_fea <- min(-1*solve(A1) %*% r_vec) > 0;
if(A0_stb == TRUE){
tot_res[[i-1]][iter, 1] <- 1;
}
if(A1_stb == TRUE){
tot_res[[i-1]][iter, 2] <- 1;
}
if(A0_fea == TRUE){
fea_res[[i-1]][iter, 1] <- 1;
}
if(A1_fea == TRUE){
fea_res[[i-1]][iter, 2] <- 1;
}
iter <- iter - 1;
}
print(i);
}
all_res <- summarise_randmat_ga(tot_res = tot_res, fea_res = fea_res);
return(all_res);
}
#' Genetic algorithm to find stability
#'
#' Runs a genetic algorithm in attempt to find a vector (gamma) that, when the
#' elements of which are multiplied as scalar values on a matrix, A1, which
#' represents a complex system, causes all real parts of eigenvalues of A1 to
#' be negative, and the complex system to therefore be stable.
#'
#'@return A value indicating whether or not a stabilising vector is (1) or is
#'not (0) found by the genetic algorithm. If so, then the vector is appended to
#'the file "evo_out.txt".
#'@param A1 Square matrix representing a complex system
#'@param max_it Maximum number of generations the genetic algorithm is allowed
#'to run before giving up and declaring no stabilising set of gammas is found
#'@param converg Convergence criteria for the genetic algorithm; if the mean
#'fitness from one generation to the next (in terms of improvement shifting
#'leading Real parts of eigenvalues to the left) is below this, then the
#'algorithm gives up and declares no stabilising set of gammas is found
#'@examples
#'A_mat <- matrix(data = rnorm(n = 16), nrow = 4);
#'sim <- rand_mat_ga(A1 = A_mat, max_it = 10);
#'@export
rand_mat_ga <- function(A1, max_it = 20, converg = 0.01){
nn <- dim(A1)[1];
rind <- runif(n = nn*1000, min = 0, max = 1);
inds <- matrix(data = rind, nrow = 1000, ncol = nn);
lastf <- -10;
ccrit <- 10;
find_st <- 0;
iter <- max_it;
while(iter > 0 & find_st < 1 & ccrit > converg){
ivar <- rep(x = 0, length = dim(inds)[1]);
ifit <- rep(x = 0, length = dim(inds)[1]);
isst <- rep(x = 0, length = dim(inds)[1]);
for(i in 1:dim(inds)[1]){
ifit[i] <- -1*max(Re(eigen(inds[i,]*A1)$values));
ivar[i] <- var(inds[i,]);
isst[i] <- max(Re(eigen(inds[i,]*A1)$values)) < 0;
}
most_fit <- order(ifit, decreasing = TRUE)[1:20];
parents <- inds[most_fit,];
new_gen <- matrix(data = t(parents), nrow = 1000, ncol = nn,
byrow = TRUE);
mu_dat <- rbinom(n = nn*1000, size = 1, prob = 0.2);
mu_dat2 <- rnorm(n = nn*1000, mean = 0, sd = 0.02);
mu_dat2[mu_dat2 < 0] <- -mu_dat2[mu_dat2 < 0];
mu_dat2[mu_dat2 > 2] <- 2;
mu_dat3 <- mu_dat * mu_dat2;
mu_mat <- matrix(data = mu_dat3, nrow = 1000, ncol = nn);
new_gen <- new_gen + mu_mat;
new_gen <- crossover(inds = new_gen, pr = 0.1);
inds <- new_gen;
find_st <- max(isst);
newf <- mean(ifit);
ccrit <- newf - lastf;
lastf <- newf;
iter <- iter - 1;
}
if(find_st == 1){
s_row <- which(isst == 1)[1];
writt <- c(nn, inds[s_row,]);
cat(writt, file = "evo_out.txt", append = TRUE);
cat("\n", file = "evo_out.txt", append = TRUE);
}
return(find_st);
}
# This function is used within the genetic algorithm
crossover <- function(inds, pr = 0.1){
crossed <- floor(dim(inds)[1] * pr);
cross1 <- sample(x = 1:dim(inds)[1], size = crossed);
cross2 <- sample(x = 1:dim(inds)[1], size = crossed);
for(i in 1:length(cross1)){
fromv <- sample(x = 1:dim(inds)[2], size = 1);
tov <- sample(x = 1:dim(inds)[2], size = 1);
temp <- inds[cross1[i],fromv:tov];
inds[cross1[i],fromv:tov] <- inds[cross2[i],fromv:tov];
inds[cross2[i],fromv:tov] <- temp;
}
return(inds);
}
#' Summarise random matrix results from genetic algorithm
#'
#' Takes the list output of the Evo_rand_gen_var function and summarises it into
#' a useable format.
#'
#'@return A table of stability and feasibility results.
#'@param tot_res The tot_res list output from the rand_gen_var function
#'@param fea_res The fea_res list output from the rand_gen_var function
#'@examples
#'sim_ga <- Evo_rand_gen_var(from = 2, to = 2, iters = 1);
#'sum_rand_ga <- summarise_randmat_ga(sim_ga$tot_res, sim_ga$fea_res);
#'@export
summarise_randmat_ga <- function(tot_res, fea_res){
sims <- length(tot_res);
all_res <- matrix(data = 0, nrow = sims, ncol = 7);
for(i in 1:sims){
A0_unst <- tot_res[[i]][,1] == FALSE;
A0_stbl <- tot_res[[i]][,1] == TRUE;
A1_unst <- tot_res[[i]][,2] == FALSE;
A1_stbl <- tot_res[[i]][,2] == TRUE;
stabled <- tot_res[[i]][,1] == FALSE & tot_res[[i]][,2] == TRUE;
unstabled <- tot_res[[i]][,1] == TRUE & tot_res[[i]][,2] == FALSE;
all_res[i, 1] <- i + 1;
all_res[i, 2] <- sum(A0_unst);
all_res[i, 3] <- sum(A0_stbl);
all_res[i, 4] <- sum(A1_unst);
all_res[i, 5] <- sum(A1_stbl);
all_res[i, 6] <- sum(stabled);
all_res[i, 7] <- sum(unstabled);
}
colnames(all_res) <- c("S", "A_unstable", "A_stable", "M_unstable",
"M_stable", "A_stabilised", "A_destabilised");
return(all_res);
}
#' Find the component response rate vectors of highest size
#'
#' Takes the output of the genetic algorithm showing the component response rate
#' (gamma) vectors found to be stabilising and returns the ones that are from
#' the largest systems (i.e., highest S)
#'
#'@return A list of vectors of component response rates (gammas)
#'@param evo_out Output of component response values
#'@param size The number of gamma vectors to be returned
#'@export
get_top_evo_out <- function(evo_out, size){
highest <- max(evo_out);
gammas <- NULL;
while(size > 0){
pos <- which(evo_out == highest);
len <- length(pos);
nli <- NULL;
for(i in 1:len){
start <- pos[i] + 1;
end <- pos[i] + highest;
gammas[[size]] <- evo_out[start:end];
size <- size - 1;
if(size == 0){
break;
}
}
highest <- highest - 1;
}
return(gammas);
}
#' Plot component response rate distributions
#'
#' Plots the component response rate distributions from the largest systems that
#' the genetic algorithm finds to be stable. At the moment, only a six panel
#' plot is allowed
#'
#'@return A list of vectors of component response rates (gammas)
#'@param evo_out Output of component response values
#'@export
plot_evo_out <- function(evo_out){
gammas <- get_top_evo_out(evo_out, size = 9);
par(mfrow = c(3, 3), mar = c(0.5, 0.5, 0.5, 0.5), oma = c(5, 5, 1, 1));
hist(gammas[[1]], main = "", breaks = 20, col = "grey", ylim = c(0, 6),
xaxt = "n", xlim = c(0, 1.1));
hist(gammas[[2]], main = "", breaks = 20, col = "grey", ylim = c(0, 6),
xaxt = "n", yaxt = "n", xlim = c(0, 1.1));
hist(gammas[[3]], main = "", breaks = 20, col = "grey", ylim = c(0, 6),
xaxt = "n", yaxt = "n", xlim = c(0, 1.1));
hist(gammas[[4]], main = "", breaks = 20, col = "grey", ylim = c(0, 6),
xaxt = "n", xlim = c(0, 1.1));
hist(gammas[[5]], main = "", breaks = 20, col = "grey", ylim = c(0, 6),
xaxt = "n", yaxt = "n", xlim = c(0, 1.1));
hist(gammas[[6]], main = "", breaks = 20, col = "grey", ylim = c(0, 6),
xaxt = "n", yaxt = "n", xlim = c(0, 1.1));
hist(gammas[[7]], main = "", breaks = 20, col = "grey", ylim = c(0, 6),
xlim = c(0, 1.1));
hist(gammas[[8]], main = "", breaks = 20, col = "grey", ylim = c(0, 6),
yaxt = "n", xlim = c(0, 1.1));
hist(gammas[[9]], main = "", breaks = 20, col = "grey", ylim = c(0, 6),
yaxt = "n", xlim = c(0, 1.1));
mtext(side = 1, outer = TRUE, line = 3, cex = 1.5,
text = expression(paste("Component response rate (",gamma,")")));
mtext(side = 2, text = "Frequency", outer = TRUE, line = 3, cex = 1.5);
}