1# These functions are
2# Copyright (C) 1998-2021 T.W. Yee, University of Auckland.
3# All rights reserved.
4
5
6
7
8
9
10
11
12fill <-
13fill1 <- fill2 <- fill3 <-
14  function(x, values = 0, ncolx = ncol(x)) {
15  x <- as.matrix(x)
16  matrix(values, nrow = nrow(x), ncol = ncolx, byrow = TRUE)
17}
18
19
20
21
22
23extract.arg <- function(a) {
24  s <- substitute(a)
25  as.character(s)
26}
27
28
29
30
31
32
33remove.arg <- function(string) {
34
35  nc <- nchar(string)
36  bits <- substring(string, 1:nc, 1:nc)
37  b1 <- (1:nc)[bits == "("]
38  b1 <- if (length(b1)) b1[1]-1 else nc
39  if (b1 == 0)
40    return("")
41  string <- paste(bits[1:b1], collapse = "")
42  string
43}
44
45
46
47
48
49
50add.arg <- function(string, arg.string) {
51
52  if (arg.string == "")
53    return(string)
54  nc <- nchar(string)
55  lastc <- substring(string, nc, nc)
56  if (lastc == ")") {
57    if (substring(string, nc-1, nc-1) == "(") {
58      paste(substring(string, 1, nc-2), "(", arg.string, ")",
59            sep = "")
60    } else {
61      paste(substring(string, 1, nc-1), ", ", arg.string, ")",
62            sep = "")
63    }
64  } else {
65    paste(string, "(", arg.string, ")", sep = "")
66  }
67}
68
69
70
71
72
73
74get.arg <- function(string) {
75
76  nc <- nchar(string)
77  bits <- substring(string, 1:nc, 1:nc)
78  b1 <- (1:nc)[bits == "("]
79  b2 <- (1:nc)[bits == ")"]
80  b1 <- if (length(b1)) min(b1) else return("")
81  b2 <- if (length(b2)) max(b2) else return("")
82  if (b2-b1 == 1) "" else paste(bits[(1+b1):(b2-1)], collapse = "")
83}
84
85
86
87
88
89
90
91 eifun <- function(i, n)
92    cbind(as.numeric((1:n) == i))
93
94
95
96 eifun <-
97 I.col <- function(i, n)
98    diag(n)[, i, drop = FALSE]
99
100
101
102 eijfun <- function(i, n) {
103  temp <- matrix(0, n, 1)
104  if (length(i))
105    temp[i, ] <- 1
106  temp
107}
108
109
110
111
112
113
114tapplymat1 <-
115  function(mat, function.arg = c("cumsum", "diff", "cumprod")) {
116
117
118  if (!missing(function.arg))
119    function.arg <- as.character(substitute(function.arg))
120  function.arg <- match.arg(function.arg,
121                            c("cumsum", "diff", "cumprod"))[1]
122
123  type <- switch(function.arg, cumsum = 1, diff = 2, cumprod = 3,
124           stop("argument 'function.arg' not matched"))
125
126  if (!is.matrix(mat))
127    mat <- as.matrix(mat)
128  NR <- nrow(mat)
129  NC <- ncol(mat)
130  fred <- .C("tapply_mat1", mat = as.double(mat), as.integer(NR),
131             as.integer(NC), as.integer(type))  # , PACKAGE = "VGAM"
132  dim(fred$mat) <- c(NR, NC)
133  dimnames(fred$mat) <- dimnames(mat)
134  switch(function.arg,
135         cumsum  = fred$mat,
136         diff    = fred$mat[, -1, drop = FALSE],
137         cumprod = fred$mat)
138}
139
140
141
142
143
144
145matrix.power <- function(wz, M, power, fast = TRUE) {
146
147
148
149
150  n <- nrow(wz)
151  index <- iam(NA, NA, M, both = TRUE, diag = TRUE)
152  dimm.value <- if (is.matrix(wz)) ncol(wz) else 1
153  if (dimm.value > M*(M+1)/2)
154    stop("too many columns")
155
156
157  if (M == 1 || dimm.value == M) {
158    WW <- wz^power  # May contain NAs
159    return(t(WW))
160  }
161
162  if (fast) {
163    k <- veigen(t(wz), M = M)  # matrix.arg)
164    evals  <- k$values   # M x n
165    evects <- k$vectors  # M x M x n
166  } else {
167    stop("sorry, cannot handle matrix-band form yet")
168    k <- unlist(apply(wz, 3, eigen), use.names = FALSE)
169    dim(k) <- c(M, M+1, n)
170    evals  <- k[,  1, , drop = TRUE]  # M x n
171    evects <- k[, -1, , drop = TRUE]  # M x M x n
172  }
173
174  temp <- evals^power    # Some values may be NAs
175
176
177  index <- as.vector( matrix(1, 1, M) %*% is.na(temp) )
178
179
180  index <- (index == 0)
181  if (!all(index)) {
182    warning("Some weight matrices have negative ",
183            "eigenvalues. They will be assigned NAs")
184    temp[,!index] <- 1
185  }
186
187  WW <- mux55(evects, temp, M = M)
188  WW[,!index] <- NA
189  WW
190}
191
192
193
194
195
196
197ResSS.vgam <- function(z, wz, M) {
198
199
200  if (M == 1)
201    return(sum(c(wz) * c(z^2)))
202
203  wz.z <- mux22(t(wz), z, M = M, as.matrix = TRUE)
204  sum(wz.z * z)
205}
206
207
208
209
210
211
212wweighted.mean <- function(y, w = NULL, matrix.arg = TRUE) {
213  if (!matrix.arg)
214    stop("currently, argument 'matrix.arg' must be TRUE")
215  y <- as.matrix(y)
216  M <- ncol(y)
217  n <- nrow(y)
218  if (M == 1) {
219    if (missing(w)) mean(y) else sum(w * y) / sum(w)
220  } else {
221    if (missing(w)) y %*% rep(1, n) else {
222      numer <- mux22(t(w), y, M, as.matrix = TRUE)
223      numer <- t(numer) %*% rep(1, n)
224      denom <- t(w) %*% rep(1, n)
225      denom <- matrix(denom, 1, length(denom))
226      if (matrix.arg)
227        denom <- m2a(denom, M = M)[, , 1]
228      c(solve(denom, numer))
229    }
230  }
231}
232
233
234
235
236
237
238veigen <- function(x, M) {
239
240
241  n <- ncol(x)
242  index <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
243  dimm.value <- nrow(x)  # usually M or M(M+1)/2
244
245  z <- .Fortran("veigenf",
246      as.integer(M),
247      as.integer(n),
248      as.double(x),
249      values = double(M * n),
250      as.integer(1),
251      vectors = double(M*M*n),
252      double(M),
253      double(M),
254      wk = double(M*M),
255      as.integer(index$row), as.integer(index$col),
256      as.integer(dimm.value),
257      error.code = integer(1))
258
259  if (z$error.code)
260    stop("eigen algorithm (rs) returned error code ", z$error.code)
261  ord <- M:1
262  dim(z$values) <- c(M, n)
263  z$values <- z$values[ord, , drop = FALSE]
264  dim(z$vectors) <- c(M, M, n)
265  z$vectors <- z$vectors[, ord, , drop = FALSE]
266  return(list(values  = z$values,
267              vectors = z$vectors))
268}
269
270
271
272
273
274
275ima <- function(j, k, M) {
276  if (length(M) > 1 || M <= 0 || j <= 0 || k <= 0 ||
277      j > M || k > M)
278    stop("input wrong in ima()")
279  m <- diag(M)
280  m[col(m) <= row(m)] <- 1:(M*(M+1)/2)
281  if (j >= k) m[j, k] else m[k, j]
282}
283
284
285
286
287
288
289checkwz <- function(wz, M, trace = FALSE,
290                    wzepsilon = .Machine$double.eps^0.75) {
291  if (wzepsilon > 0.5)
292    warning("argument 'wzepsilon' is probably too large")
293  if (!is.matrix(wz))
294    wz <- as.matrix(wz)
295  wzsubset <- wz[, 1:M, drop = FALSE]
296  if (any(is.na(wzsubset)))
297    stop("NAs found in the working weights variable 'wz'")
298  if (any(!is.finite(wzsubset)))
299    stop("Some elements in the working weights variable 'wz' are ",
300         "not finite")
301
302  if ((temp <- sum(wzsubset < wzepsilon)))
303    warning(temp, " diagonal elements of the working weights variable ",
304            "'wz' have been replaced by ", signif(wzepsilon, 5))
305  wz[, 1:M] <- pmax(wzepsilon, wzsubset)
306  wz
307}
308
309
310
311
312
313
314
315
316label.cols.y <-
317  function(answer,
318           colnames.y = NULL,
319           NOS = 1,
320           percentiles = c(25, 50, 75),
321           one.on.one = TRUE,
322           byy = TRUE) {
323  if (!is.matrix(answer))
324    answer <- as.matrix(answer)
325
326  if (one.on.one) {
327    colnames(answer) <-
328      if (length(colnames.y) == ncol(answer))
329        colnames.y else NULL
330    return(answer)
331  }
332
333
334
335  if (is.null(percentiles))
336    percentiles <- c(25, 50, 75)  # Restore to the default
337
338  if (!is.Numeric(percentiles) ||
339      min(percentiles) <= 0 ||
340      max(percentiles) >= 100)
341    stop("values of 'percentiles' should be in [0, 100]")
342
343  percentiles <- signif(percentiles, digits = 5)
344
345  ab1 <- rep(as.character(percentiles), length = ncol(answer))
346  ab1 <- paste(ab1, "%", sep = "")
347  if (NOS > 1) {
348    suffix.char <- if (length(colnames.y) == NOS)
349      colnames.y else as.character(1:NOS)
350    ab1 <- paste(ab1, rep(suffix.char, each = length(percentiles)),
351                 sep = "")
352  }
353  colnames(answer) <- ab1
354
355
356  if (byy) {
357    answer <-
358      answer[, interleave.VGAM(.M = NCOL(answer),
359                               M1 = NOS),   # length(percentiles)),
360             drop = FALSE]
361  }
362  answer
363}
364
365
366
367
368
369
370
371
372