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