/
ContactsGen-v1.1.3.R
293 lines (249 loc) · 11.8 KB
/
ContactsGen-v1.1.3.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
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
# ContactsGen-v1.1.3.R is part of Food INdustry CoViD Control Tool
# (FInd CoV Control), version 1.1.3.
# Copyright (C) 2020-2021 Cornell University.
#
# This program 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; either version 2 of the License, or
# (at your option) any later version.
#
# This program 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 this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
# Basic model: Employees are organized into crews, each with one foreman and
# one or more (ordinary) workers. Crews are organized into teams, each of which
# is under a supervisor. Supervisors, in turn, are all under a single manager.
# Crews may be of different sizes, and teams may include different numbers of crews.
# (In the current version of the complete model, all crews are the same size, and all
# teams include the same number of the crews, but this is not required by the
# functions in this file. If this requirement is not relaxed in future versions,
# this file may be simplifed by taking explicit account of it.
library(Matrix)
# The following are helper functions used in Primary_Matrices, not meant to be
# called directly by other functions
# Generates a single vector, for a single crew
Crew_Worker = function(size) {
c(0, rep(1, size)) # 0 for the foreman, 1 for the workers
}
# Generates a single vector, for a single crew
Crew_Foreman = function(size) {
c(1, rep(0, size)) # 1 for the foreman, 0 for the workers
}
# Prepends 0 to the first vector in a list, to represent the supervisor, for a
# list of crews, or the manager, for a list of teams. This is a little kludgy
# (the supervisor is not actually a member of the first crew on their team;
# the manager is not actually a member of the first team on the farm), but it
# generates the correct block matrix when fed to bdiag.
Prepend_Zero = function(list_) {
list_[[1]] = c(0, list_[[1]])
list_
}
# Generates a list of vectors for use in creating the crew_worker matrix.
# (And, after flattening into a single vector, for use in creating the
# team_worker matrix.)
# One call generates vectors for all crews in one team.
Team_Crew_Worker = function(sizes) {
Prepend_Zero(lapply(sizes, Crew_Worker)) # Prepending 0 for the supervisor
}
# Generates a list of vectors for use in creating the crew_foreman matrix.
# (And, after flattening into a single vector, for use in creating the
# team_foreman matrix.)
# One call generates vectors for all crews in one team.
Team_Crew_Foreman = function(sizes) {
Prepend_Zero(lapply(sizes, Crew_Foreman)) # Prepending 0 for the supervisor
}
# Generates a single vector, for a single team.
Team_Worker = function(sizes) {
unlist(Team_Crew_Worker(sizes), use.names = FALSE)
}
# Generates a single vector, for a single team.
Team_Foreman = function(sizes) {
unlist(Team_Crew_Foreman(sizes), use.names = FALSE)
}
# Generates a single vector, for a single team.
Team_Supervisor = function(sizes) {
c(1, rep(0, sum(sizes) + length(sizes))) # 1 for the supervisor,
# 0 for the foreman and workers.
}
#Prepends 0 to represent the manager
List_to_Matrix = function(list_) {
as.matrix(bdiag(Prepend_Zero(list_)))
}
# Generates "primary" matrices for each combination of "level" ("crew,"
# "team," or "all") and employee type (worker, foreman, supervisor,
# manager).
#
# All have N rows, one for each individual in the population.
# The number of columns varies:
# "crew" level matrices have one column per crew,
# "team" level matrices have one column per team, and
# "all" level matrices have only a single column.
#
# In all cases, a value of 1 indicates membership, and a value of 0 indicates
# lack of membership. Example: crew_worker[i,j] = 1 if and only if employee 1
# is an (ordinary) worker who is part of crew number j.
#
# For convenience, the order of employees reflects the structure of the
# workforce, i.e., employee #1 is the manager, each supervisor is followed by
# all the employees (foremen and workers) under that supervisor, and each
# foreman is followed by all the workers under that foreman.
#
# Parameters:
# crews_by_team: a vector of how many crews are in each team
# crew_sizes: a vector of how many workers (not including the foreman) are in
# each crew (ordered by team)
Primary_Matrices = function(crews_by_team, crew_sizes) {
# Turns the single vector crew_sizes into a list
# The latter could simply have been given as a parameter, rather than
# using the two current parameters. This may or may not be changed
# in the future.
crew_sizes_split = split(crew_sizes, rep(1:length(crews_by_team),
crews_by_team))
crew_worker = List_to_Matrix(unlist(lapply(crew_sizes_split,
Team_Crew_Worker),
recursive = FALSE))
crew_foreman = List_to_Matrix(unlist(lapply(crew_sizes_split,
Team_Crew_Foreman),
recursive = FALSE))
team_worker = List_to_Matrix(lapply(crew_sizes_split, Team_Worker))
team_foreman = List_to_Matrix(lapply(crew_sizes_split, Team_Foreman))
team_supervisor = List_to_Matrix(lapply(crew_sizes_split, Team_Supervisor))
all_worker = as.matrix(rowSums(team_worker))
all_foreman = as.matrix(rowSums(team_foreman))
all_supervisor = as.matrix(rowSums(team_supervisor))
# 1 for the manager, 0 for workers, foremen, and supervisors
all_manager = as.matrix(c(1, rep(0, sum(crew_sizes) + length(crew_sizes) +
length(crews_by_team))))
list(crew_worker = crew_worker,
crew_foreman = crew_foreman,
team_worker = team_worker,
team_foreman = team_foreman,
team_supervisor = team_supervisor,
all_worker = all_worker,
all_foreman = all_foreman,
all_supervisor = all_supervisor,
all_manager = all_manager)
}
# Generates "secondary" matrices for each combination of "level" ("same crew,"
# "same team [but different crew]" or "other") and two employee types (worker,
# foreman, supervisor, manager; can be the same or different).
#
# All have N rows and N columns, one for each individual in the population.
#
# For ease of future use, I've adopted a more compact notation than in
# Primary_Matrices: level is indicated by "c," "t," or "o," and employee types
# by "w," "f," "s," or "m," with an underscore separating level from types.
# Because these matrices are treated as symmetrical, only one is stored for,
# e.g., "t_wf" and "t_fw"; in all cases, the lower-level position is named first
# (e.g., "t_wf," not "t_fw").
#
# In all cases, a value of 1 indicates a pair of employees that have the
# types and share the level in question, and a value of 0 indicates a pair of
# employees that do not.
# Example: t_wf[i,j] = 1 if and only if one of employee i and employee j is a
# foreman and the other is a worker, and they are on the same team but not on
# the same crew.
Secondary_Matrices = function(primary_matrices) {
cw = primary_matrices[['crew_worker']]
tw = primary_matrices[['team_worker']]
aw = primary_matrices[['all_worker']]
cf = primary_matrices[['crew_foreman']]
tf = primary_matrices[['team_foreman']]
af = primary_matrices[['all_foreman']]
ts = primary_matrices[['team_supervisor']]
as_ = primary_matrices[['all_supervisor']] #underscore to get around builtin
am = primary_matrices[['all_manager']]
N = dim(cw)[1]
c_ww = cw %*% t(cw)
t_ww = tw %*% t(tw) - c_ww
o_ww = aw %*% t(aw) - (c_ww + t_ww)
diag(c_ww) = 0 # we do not count pairs of an individual with themself
t_ff = tf %*% t(tf)
o_ff = af %*% t(af) - t_ff
diag(t_ff) = 0 # we do not count pairs of an individual with themself
c_wf = cw %*% t(cf) + cf %*% t(cw)
t_wf = tw %*% t(tf) + tf %*% t(tw) - c_wf
o_wf = aw %*% t(af) + af %*% t(aw) - (c_wf + t_wf)
#diag inherently 0
o_ss = as_ %*% t(as_)
diag(o_ss) = 0 # we do not count pairs of an individual with themself
t_ws = tw %*% t(ts) + ts %*% t(tw)
o_ws = aw %*% t(as_) + as_ %*% t(aw) - t_ws
#diag inherently 0
t_fs = tf %*% t(ts) + ts %*% t(tf)
o_fs = af %*% t(as_) + as_ %*% t(af) - t_fs
#diag inherently 0
#o_mm is inherently 0, given that we do not count pairs of an individual
#with themself, as there is only one manager.
o_wm = aw %*% t(am) + am %*% t(aw)
o_fm = af %*% t(am) + am %*% t(af)
o_sm = as_ %*% t(am) + am %*% t(as_)
#diag inherently 0 for all three
list(c_ww = c_ww,
t_ww = t_ww,
o_ww = o_ww,
t_ff = t_ff,
o_ff = o_ff,
c_wf = c_wf,
t_wf = t_wf,
o_wf = o_wf,
o_ss = o_ss,
t_ws = t_ws,
o_ws = o_ws,
t_fs = t_fs,
o_fs = o_fs,
o_wm = o_wm,
o_fm = o_fm,
o_sm = o_sm
)
}
# Calculate work contact rates (not binary yes/no) from a list of secondary
# matrices, a list of rates for each type of pair, and an average (daily)
# contact rate _per person_ (not per pair).
# NB: As currently coded, this will give wrong results if the order of the
# matrices and the order of the rates are not the same. Care is required.
Work_Contacts = function(secondary_matrices, rates, average) {
# NB: This
raw = Reduce('+', mapply('*', secondary_matrices, rates, SIMPLIFY = FALSE))
raw * average / (mean(raw) * dim(raw)[1]) #note: could equally well use
#dim(raw)[2]; since the matrices
#are square, these are identical
}
#Since home contacts are currently being modeled as homogeneous mixing, not
#writing code to generate matrices for that.
#sample_rates here is a bunch of provisional assumptions, since we need
#*something* to use here. A future version of the model may change these
#values; it may also allow end users to specify their own.
example_rates = list(c_ww = 1, #by definition (reference value);
#these are relative rates,
#so making the rate at which a worker
#contacts another worker on the same crew 1
#may be a natural way to scale it
t_ww = 0.1,
o_ww = 0.1,
t_ff = 0.2,
o_ff = 0.1,
c_wf = 1,
t_wf = 0.1,
o_wf = 0.1,
o_ss = 0.2,
t_ws = 0.3,
o_ws = 0.1,
t_fs = 1,
o_fs = 0.1,
o_wm = 0.01,
o_fm = 0.1,
o_sm = 1
)
#Okay, now for a single master function:
ContactsGen <- function(crews_by_team, crew_sizes, contact_rates, average) {
primary_matrices = Primary_Matrices(crews_by_team, crew_sizes)
secondary_matrices = Secondary_Matrices(primary_matrices)
work_contacts = Work_Contacts(secondary_matrices, contact_rates, average)
work_contacts
}