-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsimul.R
147 lines (117 loc) · 4.76 KB
/
simul.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
library(ggplot2)
# set.seed(69)
# Main parameters
K <- c(1, 10, 100) # k values
n_1 <- n_2 <- n_3 <- 200 # initial number of workers
# Backend parameters
t_a <- 0.3 # start time
t_b <- 0.4 # end time
t_step <- 0.001 # timestep
oddness <- 30 # how radical the proposals are (i.e. the max. amount of workers that can get hired/fired per timestep)
# Misc
graph_step <- 5 # generate a health plot every nth timestep
Colors <- c("#aa55ff", "#ffaa00", "#ff007f") # colors
names(Colors) = K
health <- function(t, n) {
return(exp(-(t^4 * n - 6*t)^2) / t^4 - 5)
}
choose <- function(delta_n, delta_health, k) {
p <- runif(1,0,1)
return((delta_n >= 0 & delta_health >= 0) | (delta_n <= 0 & p < exp(delta_n/k)))
}
get_results <- function(t, n, k, death) {
delta_health <- health(t, n + proposal) - health(t, n)
if (!(choose(proposal, delta_health, k)))
proposal <- 0
n <- n + proposal
h <- health(t, n)
if (h <= 0) {
death <- m
h <- 0
} else {
plot <<- plot + geom_point(aes(n, h, col=toString(k)), size=3) +
annotate(geom="text", x=450, y=140, label=paste("step =", m))
}
return(c(proposal, h, death))
}
t <- t_a
m <- 1
death_1 <- death_2 <- death_3 <- runtime <- 2*(t_b - t_a) / t_step
stats_1 <- stats_2 <- stats_3 <- data.frame(n=rep(NA, runtime), delta_n=rep(NA, runtime), health=rep(NA, runtime), stringsAsFactors=FALSE)
while (t >= t_a) {
plot <- ggplot() + xlim(c(0,500)) + xlab("Nº de Trabalhadores") + ylim(c(0,150)) + ylab("Rentabilidade") +
stat_function(fun=health, args=list(t=t)) + labs(col="k") +
scale_color_manual(values=Colors)
proposal <- runif(1, -oddness, oddness)
if (death_1 == runtime) {
results <- get_results(t, n_1, K[1], death_1)
stats_1[m, 1] <- n_1 <- n_1 + results[1]
stats_1[m, 2] <- results[1]
stats_1[m, 3] <- results[2]
if (results[3] < death_1)
death_1 <- results[3]
} else {
stats_1[m, ] <- list(0,NA,0)
}
if (death_2 == runtime) {
results <- get_results(t, n_2, K[2], death_2)
stats_2[m, 1] <- n_2 <- n_2 + results[1]
stats_2[m, 2] <- results[1]
stats_2[m, 3] <- results[2]
if (results[3] < death_2)
death_2 <- results[3]
} else {
stats_2[m, ] <- list(0,NA,0)
}
if (death_3 == runtime) {
results <- get_results(t, n_3, K[3], death_3)
stats_3[m, 1] <- n_3 <- n_3 + results[1]
stats_3[m, 2] <- results[1]
stats_3[m, 3] <- results[2]
if (results[3] < death_3)
death_3 <- results[3]
} else {
stats_3[m, ] <- list(0,NA,0)
}
# if (all(c(death_1, death_2, death_3) != runtime))
# stop("Game over")
if ((m - 1) %% graph_step == 0)
ggsave(paste("plot", m, ".jpg", sep=""))
m <- m + 1
t <- t + t_step
if (isTRUE(all.equal(t, t_b)))
t_step <- -t_step
}
ggplot() + xlab("Fluxo de Trabalhadores") + ylab("Nº de Ocorrências") + xlim(c(-oddness, oddness)) +
geom_histogram(aes(x=stats_1$delta_n), binwidth=1)
ggsave("fluxo_1.jpg")
ggplot() + xlab("Fluxo de Trabalhadores") + ylab("Nº de Ocorrências") + xlim(c(-oddness, oddness)) +
geom_histogram(aes(x=stats_2$delta_n), binwidth=1)
ggsave("fluxo_2.jpg")
ggplot() + xlab("Fluxo de Trabalhadores") + ylab("Nº de Ocorrências") + xlim(c(-oddness, oddness)) +
geom_histogram(aes(x=stats_3$delta_n), binwidth=1)
ggsave("fluxo_3.jpg")
ggplot() + labs(col = "k") + xlab("Tempo") + ylab("Emprego Marginal") +
xlim(c(0,min(death_1, death_2, death_3))) + ylim(c(-oddness, oddness)) +
geom_line(aes(x=as.numeric(row.names(stats_1)), y=stats_1$delta_n, col=toString(K[1]))) +
geom_line(aes(x=as.numeric(row.names(stats_2)), y=stats_2$delta_n, col=toString(K[2]))) +
geom_line(aes(x=as.numeric(row.names(stats_3)), y=stats_3$delta_n, col=toString(K[3]))) +
scale_color_manual(values=Colors)
ggsave("emprego.jpg")
ggplot() + labs(col = "k") + xlab("Tempo") + ylab("Nº de Trabalhadores") +
geom_line(data=stats_1, aes(x=as.numeric(row.names(stats_1)), y=n, col=toString(K[1]))) +
geom_line(data=stats_2, aes(x=as.numeric(row.names(stats_2)), y=n, col=toString(K[2]))) +
geom_line(data=stats_3, aes(x=as.numeric(row.names(stats_3)), y=n, col=toString(K[3]))) + scale_color_manual(values=Colors)
ggsave("trabalhadores.jpg")
ggplot() + labs(col = "k") + xlab("Tempo") + ylab("Rentabilidade") +
geom_line(data=stats_1, aes(x=as.numeric(row.names(stats_1)), y=health, col=toString(K[1]))) +
geom_line(data=stats_2, aes(x=as.numeric(row.names(stats_2)), y=health, col=toString(K[2]))) +
geom_line(data=stats_3, aes(x=as.numeric(row.names(stats_3)), y=health, col=toString(K[3]))) + scale_color_manual(values=Colors)
ggsave("rentabilidade.jpg")
d <- data.frame(k=paste(unlist(K)), lifespan=c(death_1, death_2, death_3))
ggplot(d, aes(x=k, y=lifespan, label=lifespan, fill=k)) + xlab("k") + ylab("Iterações Sobrevividas") +
geom_bar(stat="identity") +
geom_text(size = 9) +
scale_fill_manual(values=Colors) +
scale_x_discrete(limits=paste(unlist(K)))
ggsave("vida.jpg")