1# These functions are 2# Copyright (C) 1998-2021 T.W. Yee, University of Auckland. 3# All rights reserved. 4 5 6 7 8mux34 <- function(xmat, cc, symmetric = FALSE) { 9 10 11 if (!is.matrix(xmat)) 12 xmat <- as.matrix(xmat) 13 d <- dim(xmat) 14 nnn <- d[1] 15 RRR <- d[2] 16 if (length(cc) == 1) 17 cc <- matrix(cc, 1, 1) 18 if (!is.matrix(cc)) 19 stop("'cc' is not a matrix") 20 c( .C("VGAM_C_mux34", as.double(xmat), as.double(cc), 21 as.integer(nnn), as.integer(RRR), 22 as.integer(symmetric), ans = as.double(rep_len(0.0, nnn)), 23 NAOK = TRUE)$ans) 24} 25 26 27 28 29 30 31 32 33 34mux2 <- function(cc, xmat) { 35 36 37 if (!is.matrix(xmat)) 38 xmat <- as.matrix(xmat) 39 d <- dim(xmat) 40 n <- d[1] 41 p <- d[2] 42 if (is.matrix(cc)) 43 cc <- array(cc, c(dim(cc), n)) 44 d <- dim(cc) 45 M <- d[1] 46 if (d[2] != p || d[3] != n) 47 stop("dimension size inconformable") 48 ans <- rep_len(NA_real_, n*M) 49 fred <- .C("mux2ccc", as.double(cc), as.double(t(xmat)), 50 ans = as.double(ans), as.integer(p), as.integer(n), 51 as.integer(M), NAOK = TRUE) 52 matrix(fred$ans, n, M, byrow = TRUE) 53} 54 55 56 57 58 59 60mux22 <- function(cc, xmat, M, upper = FALSE, as.matrix = FALSE) { 61 62 n <- ncol(cc) 63 64 index <- iam(NA, NA, M, both = TRUE, diag = TRUE) 65 dimm.value <- nrow(cc) # Usually M or M(M+1)/2 66 67 ans <- rep_len(NA_real_, n*M) 68 fred <- .C("mux22ccc", as.double(cc), as.double(t(xmat)), 69 ans = as.double(ans), as.integer(dimm.value), 70 as.integer(index$row), as.integer(index$col), 71 as.integer(n), as.integer(M), wk = double(M*M), 72 as.integer(as.numeric(upper)), NAOK = TRUE) 73 if (as.matrix) { 74 dim(fred$ans) <- c(M, n) 75 t(fred$ans) 76 } else { 77 fred$ans 78 } 79} 80 81 82 83 84 85mux5 <- function(cc, x, M, matrix.arg = FALSE) { 86 87 88 89 dimx <- dim(x) 90 dimcc <- dim(cc) 91 r <- dimx[2] 92 93 if (matrix.arg) { 94 n <- dimcc[1] 95 neltscci <- ncol(cc) 96 cc <- t(cc) 97 } else { 98 n <- dimcc[3] 99 if (dimcc[1] != dimcc[2] || 100 dimx[1] != dimcc[1] || 101 (length(dimx) == 3 && dimx[3] != dimcc[3])) 102 stop('input nonconformable') 103 neltscci <- M*(M+1)/2 104 } 105 106 if (is.matrix(x)) 107 x <- array(x, c(M, r, n)) 108 index.M <- iam(NA, NA, M, both = TRUE, diag = TRUE) 109 index.r <- iam(NA, NA, r, both = TRUE, diag = TRUE) 110 111 size <- if (matrix.arg) dimm(r) * n else r * r * n 112 fred <- .C("mux5ccc", as.double(cc), as.double(x), 113 ans = double(size), 114 as.integer(M), as.integer(n), as.integer(r), 115 as.integer(neltscci), 116 as.integer(dimm(r)), 117 as.integer(as.numeric(matrix.arg)), 118 double(M*M), double(r*r), 119 as.integer(index.M$row), as.integer(index.M$col), 120 as.integer(index.r$row), as.integer(index.r$col), 121 ok3 = as.integer(1), NAOK = TRUE) 122 if (fred$ok3 == 0) 123 stop("can only handle 'matrix.arg == 1'") 124 125 126 if (matrix.arg) { 127 ans <- fred$ans 128 dim(ans) <- c(dimm(r), n) 129 t(ans) 130 } else { 131 array(fred$ans, c(r, r, n)) 132 } 133} 134 135 136 137 138mux55 <- function(evects, evals, M) { 139 140 d <- dim(evects) 141 n <- ncol(evals) 142 if (d[1] != M || d[2] != M || d[3] != n || 143 nrow(evals)!= M || ncol(evals) != n) 144 stop("input nonconformable") 145 MMp1d2 <- M*(M+1)/2 # The answer is a full-matrix 146 index <- iam(NA, NA, M, both = TRUE, diag = TRUE) 147 148 fred <- .C("mux55ccc", as.double(evects), as.double(evals), 149 ans = double(MMp1d2 * n), 150 double(M*M), double(M*M), 151 as.integer(index$row), as.integer(index$col), 152 as.integer(M), as.integer(n), NAOK = TRUE) 153 dim(fred$ans) <- c(MMp1d2, n) 154 fred$ans 155} 156 157 158 159 160mux7 <- function(cc, x) { 161 162 dimx <- dim(x) 163 dimcc <- dim(cc) 164 if (dimx[1]!= dimcc[2] || 165 (length(dimx) == 3 && dimx[3]!= dimcc[3])) 166 stop('input nonconformable') 167 M <- dimcc[1] 168 qq <- dimcc[2] 169 n <- dimcc[3] 170 r <- dimx[2] 171 if (is.matrix(x)) 172 x <- array(x, c(qq, r, n)) 173 174 ans <- array(NA, c(M, r, n)) 175 fred <- .C("mux7ccc", as.double(cc), as.double(x), 176 ans = as.double(ans), 177 as.integer(M), as.integer(qq), as.integer(n), 178 as.integer(r), NAOK = TRUE) 179 array(fred$ans, c(M, r, n)) 180} 181 182 183 184 185 186 187 188 189 190mux11 <- function(cc, xmat) { 191 192 193 dcc <- dim(cc) 194 d <- dim(xmat) 195 M <- dcc[1] 196 R <- d[2] 197 n <- dcc[3] 198 if (M != dcc[2] || d[1] != n*M) 199 stop("input inconformable") 200 201 Xmat <- array(c(t(xmat)), c(R, M, n)) 202 Xmat <- aperm(Xmat, c(2, 1, 3)) # Xmat becomes M x R x n 203 mat <- mux7(cc, Xmat) # mat is M x R x n 204 mat <- aperm(mat, c(2, 1, 3)) # mat becomes R x M x n 205 mat <- matrix(c(mat), n*M, R, byrow = TRUE) 206 mat 207} 208 209 210 211mux111 <- function(cc, xmat, M, upper = TRUE) { 212 213 214 if (!is.matrix(xmat)) 215 xmat <- as.matrix(xmat) 216 if (!is.matrix(cc)) 217 cc <- t(as.matrix(cc)) 218 R <- ncol(xmat) 219 n <- nrow(xmat) / M 220 index <- iam(NA, NA, M, both = TRUE, diag = TRUE) 221 dimm.value <- nrow(cc) # Usually M or M(M+1)/2 222 223 fred <- .C("mux111ccc", 224 as.double(cc), b = as.double(t(xmat)), 225 as.integer(M), as.integer(R), as.integer(n), 226 wkcc = double(M * M), wk2 = double(M * R), 227 as.integer(index$row), as.integer(index$col), 228 as.integer(dimm.value), 229 as.integer(upper), NAOK = TRUE) 230 231 ans <- fred$b 232 dim(ans) <- c(R, n * M) 233 d <- dimnames(xmat) 234 dimnames(ans) <- list(d[[2]], d[[1]]) 235 t(ans) 236} 237 238 239 240 241 242 243mux15 <- function(cc, xmat) { 244 245 n <- nrow(xmat) 246 M <- ncol(xmat) 247 if (nrow(cc) != M || ncol(cc) != M) 248 stop("input inconformable") 249 if (max(abs(t(cc)-cc))>0.000001) 250 stop("argument 'cc' is not symmetric") 251 252 ans <- rep_len(NA_real_, n*M*M) 253 fred <- .C("mux15ccc", as.double(cc), as.double(t(xmat)), 254 ans = as.double(ans), as.integer(M), 255 as.integer(n), NAOK = TRUE) 256 array(fred$ans, c(M, M, n)) 257} 258 259 260 261 262 263vforsub <- function(cc, b, M, n) { 264 265 266 267 index <- iam(NA, NA, M, both = TRUE, diag = TRUE) 268 dimm.value <- nrow(cc) # M or M(M+1)/2 269 270 271 fred <- .C("vforsubccc", as.double(cc), b = as.double(t(b)), 272 as.integer(M), as.integer(n), wk = double(M*M), 273 as.integer(index$row), as.integer(index$col), 274 as.integer(dimm.value), NAOK = TRUE) 275 276 dim(fred$b) <- c(M, n) 277 fred$b 278} 279 280 281 282 283vbacksub <- function(cc, b, M, n) { 284 index <- iam(NA, NA, M, both = TRUE, diag = TRUE) 285 dimm.value <- nrow(cc) 286 if (nrow(b) != M || ncol(b) != n) 287 stop("dimension size inconformable") 288 289 fred <- .C("vbacksubccc", as.double(cc), b = as.double(b), 290 as.integer(M), as.integer(n), wk = double(M*M), 291 as.integer(index$row), as.integer(index$col), 292 as.integer(dimm.value), NAOK = TRUE) 293 294 if (M == 1) { 295 fred$b 296 } else { 297 dim(fred$b) <- c(M, n) 298 t(fred$b) 299 } 300} 301 302 303 304vchol <- function(cc, M, n, silent = FALSE, callno = 0) { 305 306 307 308 309 310 index <- iam(NA, NA, M = M, both = TRUE, diag = TRUE) 311 cc <- t(cc) 312 MM <- nrow(cc) # cc is big enough to hold its Cholesky decom. 313 314 fred <- .C("vcholccc", cc = as.double(cc), as.integer(M), 315 as.integer(n), ok = integer(n), 316 wk = double(M*M), as.integer(index$row), 317 as.integer(index$col), 318 as.integer(MM), 319 NAOK = TRUE) 320 321 failed <- (fred$ok != 1) 322 if ((correction.needed <- any(failed))) { 323 index <- (1:n)[failed] 324 if (!silent) { 325 if (length(index) < 11) 326 warning("weight matri", 327 ifelse(length(index) > 1, "ces ","x "), 328 paste(index, collapse = ", "), 329 " not positive-definite") 330 } 331 } 332 333 ans <- fred$cc 334 dim(ans) <- c(MM, n) 335 336 if (correction.needed) { 337 temp <- cc[, index, drop = FALSE] 338 tmp777 <- vchol.greenstadt(temp, M = M, silent = silent, 339 callno = callno + 1) 340 341 342 if (length(index) == n) { 343 ans <- tmp777[1:nrow(ans), , drop = FALSE] 344 } else { 345 346 347 ans[, index] <- tmp777 # restored 20031016 348 } 349 } 350 dim(ans) <- c(MM, n) # Make sure 351 352 ans 353} 354 355 356 357vchol.greenstadt <- function(cc, M, silent = FALSE, 358 callno = 0) { 359 360 361 362 363 MM <- dim(cc)[1] 364 n <- dim(cc)[2] 365 366 if (!silent) 367 cat(paste("Applying Greenstadt modification to ", n, " matri", 368 ifelse(n > 1, "ces", "x"), "\n", sep = "")) 369 370 371 372 373 374 temp <- veigen(cc, M = M) # , mat = TRUE) 375 dim(temp$vectors) <- c(M, M, n) # Make sure (when M = 1) for mux5 376 dim(temp$values) <- c(M, n) # Make sure (when M = 1) for mux5 377 378 is.neg <- (temp$values < .Machine$double.eps) 379 is.pos <- (temp$values > .Machine$double.eps) 380 zilch <- (!is.pos & !is.neg) 381 382 temp$values <- abs(temp$values) 383 384 temp.small.value <- quantile(temp$values[!zilch], prob = 0.15) 385 if (callno > 2) { 386 temp.small.value <- abs(temp.small.value) * 1.50^callno 387 388 389 small.value <- temp.small.value 390 391 392 temp$values[zilch] <- small.value 393 } 394 395 if (callno > 9) { 396 warning("taking drastic action; setting all wz to ", 397 "scaled versions of the order-M identity matrix") 398 399 cc2mean <- abs(colMeans(cc[1:M, , drop = FALSE])) 400 temp$values <- matrix(cc2mean, M, n, byrow = TRUE) 401 temp$vectors <- array(c(diag(M)), c(M, M, n)) 402 } 403 404 405 406 temp3 <- mux55(temp$vectors, temp$values, M = M) #, matrix.arg = TRUE) 407 ans <- vchol(t(temp3), M = M, n = n, silent = silent, 408 callno = callno + 1) #, matrix.arg = TRUE) 409 410 411 412 if (nrow(ans) == MM) ans else ans[1:MM, , drop = FALSE] 413} 414 415 416 417 418 419 420 421 422