/
misc.R
428 lines (415 loc) · 13.9 KB
/
misc.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
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
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
# Copyright (c) 2022 Merck & Co., Inc., Rahway, NJ, USA and its affiliates. All rights reserved.
#
# This file is part of the gMCPLite program.
#
# gMCPLite 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 3 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, see <http://www.gnu.org/licenses/>.
checkValidWeights <- function(weights) {
if(!is.numeric(weights)) {
stop("Weights have to be numeric!")
}
if(any(is.na(weights) | is.infinite(weights))) {
warning("Some of the weights are not real numbers. NA, NaN, Inf and -Inf are not supported.")
}
if(any(0 > weights | weights > 1 + .Machine$double.eps ^ 0.25)) {
warning("Invalid weights: weights must be between 0 and 1")
}
if(sum(weights) > 1 + .Machine$double.eps ^ 0.25) {
warning("Invalid weights: the sum of all weights must be less than 1")
}
}
#' Substitute Epsilon
#'
#' Substitute Epsilon with a given value.
#'
#' For details see the given references.
#'
#' @param graph A graph of class \code{graphMCP} or class
#' \code{entangledMCP}.
#' @param eps A numeric scalar specifying a value for epsilon edges.
#' @return A graph where all epsilons have been replaced with the given value.
#' @author Kornelius Rohmeyer \email{rohmeyer@@small-projects.de}
#' @seealso \code{graphMCP}, \code{entangledMCP}
#' @keywords print graphs
#' @examples
#'
#'
#' graph <- improvedParallelGatekeeping()
#' graph
#' substituteEps(graph, eps=0.01)
#'
#'
#' @export substituteEps
substituteEps <- function(graph, eps=10^(-3)) {
# Call this function recursivly for entangled graphs.
if ("entangledMCP" %in% class(graph)) {
for(i in 1:length(graph@subgraphs)) {
graph@subgraphs[[i]] <- substituteEps(graph@subgraphs[[i]], eps)
}
return(graph)
}
# Real function:
if (is.numeric(graph@m)) return(graph)
m <- matrix(gsub("\\\\epsilon", eps, graph@m), nrow=length(getNodes(graph)))
m2 <- matrix(sapply(m, function(x) {
result <- try(eval(parse(text=x)), silent=TRUE);
ifelse(class(result)=="try-error",NA,result)
}), nrow=length(getNodes(graph)))
if (all(is.na(m)==is.na(m2))) m <- m2
rownames(m) <- colnames(m) <- getNodes(graph)
graph@m <- m
return(graph)
}
#' Replaces variables in a general graph with specified numeric values
#'
#' Given a list of variables and real values a general graph is processed and
#' each variable replaced with the specified numeric value.
#'
#'
#' @param graph A graph of class \code{graphMCP} or class
#' \code{entangledMCP}.
#' @param variables A named list with one or more specified real values, for example
#' \code{list(a=0.5, b=0.8, "tau"=0.5)} or \code{list(a=c(0.5, 0.8), b=0.8, "tau"=0.5)}.
#' If \code{ask=TRUE} and this list is
#' missing at all or single variables are missing from the list, the user is
#' asked for the values (if the session is not interactive an error is thrown).
#' For interactively entered values only single numbers are supported.
#' @param ask If \code{FALSE} all variables that are not specified are not
#' replaced.
#' @param partial IF \code{TRUE} only specified variables are replaced and
#' parameter \code{ask} is ignored.
#' @param expand Used internally. Don't use yourself.
#' @param list If \code{TRUE} the result will always be a list, even if only one
#' graph is returned in this list.
#' @return A graph or a matrix with variables replaced by the specified numeric
#' values. Or a list of theses graphs and matrices if a variable had more than one value.
#' @author Kornelius Rohmeyer \email{rohmeyer@@small-projects.de}
#' @seealso \code{graphMCP}, \code{entangledMCP}
#' @keywords print graphs
#' @examples
#'
#' graph <- HungEtWang2010()
#' replaceVariables(graph, list("tau"=0.5,"omega"=0.5, "nu"=0.5))
#' replaceVariables(graph, list("tau"=c(0.1, 0.5, 0.9),"omega"=c(0.2, 0.8), "nu"=0.4))
#'
#' @export replaceVariables
replaceVariables <-function(graph, variables=list(), ask=TRUE, partial=FALSE, expand=TRUE, list=FALSE) {
if (expand) variables <- varcombs(variables)
if (is.list(variables[[1]])) {
result <- list()
for (v in variables) {
r2 <- replaceVariables(graph, variables=v, ask=ask, partial=partial, expand=FALSE)
attr(r2, "label") <- attr(v, "label")
attr(r2, "variables") <- unlist(v)
result <- c(result, list(r2))
}
if (length(result)==1 && !list) return(result[[1]])
return(result)
}
# Call this function recursivly for entangled graphs.
if ("entangledMCP" %in% class(graph)) {
for(i in 1:length(graph@subgraphs)) {
graph@subgraphs[[i]] <- replaceVariables(graph@subgraphs[[i]], variables=variables, ask=ask, partial=partial)
}
return(graph)
}
# Real function:
greek <- c("alpha", "beta", "gamma", "delta", "epsilon", "zeta", "eta",
"theta", "iota", "kappa", "lambda", "mu", "nu", "xi",
"omicron", "pi", "rho", "sigma", "tau", "nu", "phi",
"chi", "psi", "omega")
if (is.matrix(graph)) { m <- graph } else {m <- graph@m}
for (g in c(greek, letters)) {
if (length(grep(g, m))!=0) {
if (is.null(answer <- variables[[g]])) {
if (!partial && ask) {
if(interactive()) {
answer <- readline(paste("Value for variable ",g,"? ", sep=""))
} else {
stop(paste("Value for variable",g,"not specified."))
}
}
}
if(!is.null(answer)) {
m <- gsub(paste(ifelse(nchar(g)==1,"","\\\\"), g, sep=""), answer, m)
}
}
}
if (is.matrix(graph)) return(parse2numeric(m))
graph@m <- m
return(parse2numeric(graph))
}
# Parses matrices of graphs (simple and entangled)
#
# Parses matrices of graphs (simple and entangled) when values are of type character, e.g. "1/3".
#
# @param g Graph of class \code{graphMCP} or \code{entangledMCP}.
# @param force Logical whether conversion to numeric should be forced or not.
# If forced all values that could not be parsed will be \code{NA}.
# Otherwise the original unchanged graph will be returned.
# @author Kornelius Rohmeyer \email{rohmeyer@@small-projects.de}
# @keywords Converted graph (if all values could be parsed or \code{force=TRUE}) or original graph.
# @examples
#
# # Nothing changes:
# gMCP:::parse2numeric(HungEtWang2010())
# # Note that other methods like printing don't handle NAs well:
# gMCP:::parse2numeric(HungEtWang2010(), force=TRUE)
#
parse2numeric <- function(graph, force=FALSE) {
# Call this function recursivly for entangled graphs.
if ("entangledMCP" %in% class(graph)) {
for(i in 1:length(graph@subgraphs)) {
graph@subgraphs[[i]] <- parse2numeric(graph@subgraphs[[i]])
}
return(graph)
}
# Real function:
if (is.matrix(graph)) { m <- graph } else {m <- graph@m}
names <- rownames(m)
m <- matrix(sapply(m, function(x) {
result <- try(eval(parse(text=x)), silent=TRUE);
ifelse(class(result)=="try-error",NA,result)
}), nrow=dim(m)[1])
if (!force && any(is.na(m))) return(graph)
rownames(m) <- colnames(m) <- names
if (is.matrix(graph)) return(m)
graph@m <- m
return(graph)
}
isEpsilon <- function(w) {
x <- try(eval(parse(text = gsub("\\\\epsilon", 0, w)), envir = baseenv()), silent=TRUE)
if ("try-error" %in% class(x)) return(FALSE)
return(x==0)
}
# For testing purposes: variables <- list(a=c(1,2), b=(3), x=c(2,3,4), d=c(1,2))
varcombs <- function(variables) {
combs <- list()
m <- do.call(expand.grid, lapply(variables, function(x){1:length(x)}))
for (i in 1:dim(m)[1]) {
variablesII <- rep(0, length(variables))
for(k in 1:length(variables)) {
variablesII[k] <- variables[[k]][m[i,k]]
}
names(variablesII) <- names(variables)
x <- as.list(variablesII)
attr(x, "label") <- paste(paste(names(variables),"=",variablesII,sep=""), collapse=", ")
combs <- c(combs, list(x))
}
# GII <- replaceVariables(G, as.list(variablesII))
# additionalLabel <- paste(",", paste(paste(names(variables),"=",variablesII,sep=""), collapse=", "))
return(combs)
}
#' Create a Block Diagonal Matrix with NA outside the diagonal
#'
#' Build a block diagonal matrix with NA values outside the diagonal given
#' several building block matrices.
#'
#' This function is useful to build the correlation matrices, when only
#' partial knowledge of the correlation exists.
#'
#' @param ... individual matrices or a \code{list} of matrices.
#' @return A block diagonal matrix with NA values outside the diagonal.
#' @author Kornelius Rohmeyer \email{rohmeyer@@small-projects.de}
#' @seealso \code{\link{gMCP}}
#' @examples
#' bdiagNA(diag(3), matrix(1/2,nr=3,nc=3), diag(2))
#'
#' @export bdiagNA
bdiagNA <- function(...) {
if (nargs() == 0)
return(matrix(nrow=0, ncol=0))
if (nargs() == 1 && !is.list(...))
return(as.matrix(...))
asList <- if (nargs() == 1 && is.list(...)) list(...)[[1]] else list(...)
if (length(asList) == 1)
return(as.matrix(asList[[1]]))
n <- 0
for (m in asList) {
if (!is.matrix(m)) {
stop("Only matrices are allowed as arguments.")
}
if (dim(m)[1]!=dim(m)[2]) {
stop("Only quadratic matrices are allowed.")
}
n <- n + dim(m)[1]
}
M <- matrix(NA, nrow=n, ncol=n)
k <- 0
for (m in asList) {
for (i in 1:dim(m)[1]) {
for (j in 1:dim(m)[1]) {
M[i+k,j+k] <- m[i,j]
}
}
k <- k + dim(m)[1]
}
return(M)
}
#' Check correlation matrix
#'
#' Sanity checks for the correlation matrix.
#'
#' @param m TBA
#' @param returnMessage TBA
#' @param na.allowed TBA
#'
#' @details
#' Checks the following properties:
#' \itemize{
#' \item Values must be between -1 and 1.
#' \item Diagonal must be equal to 1.
#' \item Matrix must be symmetric.
#' }
#'
#' @return Logical
#'
#' @export checkCorrelation
#'
#' @examples
#' NULL
checkCorrelation <- function(m, returnMessage = FALSE, na.allowed = TRUE) {
if (!na.allowed && any(is.na(m))) {
if (returnMessage) {
return("Matrix can not contain NAs.")
}
return(FALSE)
}
if (!is.numeric(m) || !is.matrix(m)) {
if (returnMessage) {
return("Matrix must be a numeric matirx.")
}
return(FALSE)
}
if (!isTRUE(all.equal(1, max(1, max(abs(m)[!is.na(m)]))))) {
if (returnMessage) {
return("Values must be between -1 and 1.")
}
return(FALSE)
}
if (!isTRUE(all.equal(diag(m), rep(1, dim(m)[1]), check.attributes = FALSE))) {
if (returnMessage) {
return("Diagonal must be equal to 1.")
}
return(FALSE)
}
if (!isSymmetric(unname(m))) {
if (returnMessage) {
return("Matrix must be symmetric.")
}
return(FALSE)
}
return(TRUE)
}
# Calculation time note: n=22 needs 12 seconds on my computer.
# With each further step calculation time nearly doubles.
#' Permutation for a design matrix
#'
#' @param n dimension of the matrix
#'
#' @return a n*(2^n) dimensional matrix
#'
#' @examples
#' permutations(3)
#'
#' @export permutations
#'
permutations <- function(n) {
outer(
(1:(2^n)) - 1, (n:1) - 1,
FUN = function(x, y) {
(x %/% 2^y) %% 2
}
)
}
#' Placement of graph nodes
#'
#' Places the nodes of a graph according to a specified layout.
#'
#' If one of \code{nrow} or \code{ncol} is not given, an attempt is made to
#' infer it from the number of nodes of the \code{graph} and the other
#' parameter. If neither is given, the graph is placed as a circle.
#'
#' @param graph A graph of class \code{graphMCP} or class
#' \code{entangledMCP}.
#' @param nrow The desired number of rows.
#' @param ncol The desired number of columns.
#' @param byrow Logical whether the graph is filled by rows (otherwise by
#' columns).
#' @param topdown Logical whether the rows are filled top-down or bottom-up.
#' @param force Logical whether a graph that has already a layout should be
#' given the specified new layout.
#' @return The graph with nodes placed according to the specified layout.
#' @author Kornelius Rohmeyer \email{rohmeyer@@small-projects.de}
#' @seealso \code{graphMCP}, \code{entangledMCP}
#' @keywords graphs
#' @examples
#'
#'
#' g <- matrix2graph(matrix(0, nrow=6, ncol=6))
#'
#' g <- placeNodes(g, nrow=2, force=TRUE)
#'
#'
#' @export placeNodes
placeNodes <- function(graph, nrow, ncol, byrow = TRUE, topdown = TRUE, force = FALSE) {
entangled <- NULL
# If the graph is entangled only place the nodes of the first graph
if ("entangledMCP" %in% class(graph)) {
entangled <- graph
graph <- entangled@subgraphs[[1]]
}
# Only place nodes if no placement data exists or parameter force is set to TRUE
if (is.null(graph@nodeAttr$X) || force) {
n <- length(getNodes(graph))
if (missing(nrow) && missing(ncol)) {
v <- (1:n)/n*2*pi
nodeX <- 300 + 250*sin(v)
nodeY <- 300 + 250*cos(v)
} else {
if (missing(nrow)) {
nrow <- ceiling(length(getNodes(graph))/ncol)
}
if (missing(ncol)) {
ncol <- ceiling(length(getNodes(graph))/nrow)
}
# if nrow*ncol<n increase the number of rows
if (nrow*ncol<n) {
nrow <- ceiling(length(getNodes(graph))/ncol)
}
if (byrow) {
nodeX <- rep(((1:ncol)-1)*200+100, nrow)
nodeY <- rep(((1:nrow)-1)*200+100, each = ncol)
} else {
nodeX <- rep(((1:ncol)-1)*200+100, each = nrow)
nodeY <- rep(((1:nrow)-1)*200+100, ncol)
}
if (!topdown) nodeY <- max(nodeY) - nodeY + 100
}
graph@nodeAttr$X <- nodeX[1:n]
graph@nodeAttr$Y <- nodeY[1:n]
for (i in getNodes(graph)) {
for (j in getNodes(graph)) {
if (graph@m[i,j]!=0) {
edgeAttr(graph, i, j, "labelX") <- -100
edgeAttr(graph, i, j, "labelY") <- -100
}
}
}
}
if (!is.null(entangled)) {
entangled@subgraphs[[1]] <- graph
return(entangled)
}
return(graph)
}