1# These functions are
2# Copyright (C) 1998-2021 T.W. Yee, University of Auckland.
3# All rights reserved.
4
5
6
7
8 psv2magic <-
9    function(x.VLM, constraints, spar.vlm, sm.osps.list) {
10
11
12
13
14  colperm <- function(x, from, to) {
15
16    ncx <- ncol(x)
17    if (length(from) != length(to) ||
18        any(from != round(from)) ||
19        any(from < 1 | ncx < from) ||
20        any(duplicated(from)) ||
21        any(sort(from) != sort(to)))
22      stop("invalid column permutation indices")
23    perm <- seq_len(ncx)
24    perm[to] <- perm[from]
25    x[, perm]
26  }
27
28
29
30  assignx <- sm.osps.list$assignx
31  nassignx <- names(assignx)
32  indexterms <- sm.osps.list$indexterms
33  which.X.sm.osps <- sm.osps.list$which.X.sm.osps
34  term.labels <- sm.osps.list$term.labels
35  ncol.X.sm.osps <- sapply(which.X.sm.osps, length)
36  ncolHlist.model <- unlist(lapply(constraints, ncol))
37
38
39  ncolHlist.new <- ncolHlist.model
40  if (names(constraints)[[1]] == "(Intercept)") {
41    ncolHlist.new <- ncolHlist.new[-1]
42    nassignx <- nassignx[-1]
43  }
44
45
46  ncol.H.ps <- ncolHlist.new[indexterms]
47  num.osps.terms <- length(which.X.sm.osps)
48
49
50  ncol.allterms <- sapply(assignx, length)
51
52  ncol.model <- if (names(constraints)[[1]] == "(Intercept)")
53                ncol.allterms[-1] else ncol.allterms
54  jay <- 0
55  jjoffset <- if (names(constraints)[[1]] == "(Intercept)")
56              ncolHlist.model[1] else 0
57  perm.list <- list()
58  for (ii in seq_along(term.labels)) {
59    if (indexterms[ii]) {
60      jay <- jay + 1
61      perm.list[[jay]] <-
62        matrix(jjoffset + 1:(ncol.X.sm.osps[jay] * ncol.H.ps[jay]),
63               nrow = ncol.X.sm.osps[jay],  # Redundant really
64               ncol = ncol.H.ps[jay], byrow = TRUE)
65      jjoffset <- jjoffset + ncol.H.ps[[jay]] * ncol.X.sm.osps[[jay]]
66    } else {
67      jjoffset <- jjoffset + ncolHlist.new[ii] * ncol.model[ii]
68    }
69  }  # for ii
70  vindex.min <- sapply(perm.list, min)  # function(x) min(x)
71  vindex.max <- sapply(perm.list, max)  # function(x) max(x)
72  oo1 <- vector("list", length(ncol.H.ps))  # list()
73  for (ii in seq_along(ncol.H.ps)) {
74    oo1[[ii]] <- seq.int(vindex.min[ii], vindex.max[ii])
75  }
76  ooo <- unlist(oo1, use.names = FALSE)  # do.call("c", oo1)
77  ppp <- unlist(perm.list, use.names = FALSE)  # do.call("c", perm.list)
78
79
80  OFF.list <- vector("list", num.osps.terms)  # list()
81  for (ii in 1:num.osps.terms) {
82    index <- 0
83    OFF.list[[ii]] <- numeric()
84    for (jay in 1:(ncol.H.ps[ii])) {
85      OFF.list[[ii]][jay] <- vindex.min[ii] + index
86      index <- ncol.X.sm.osps[ii] * jay
87    }
88  }
89
90
91  list(x.VLM.new = if (identical(ppp, ooo)) x.VLM else
92                   colperm(x.VLM, ppp, ooo),
93       sp = unlist(spar.vlm),
94       S.arg = rep(sm.osps.list$S.arg, ncol.H.ps),  # Argument 'S' of magic()
95       OFF = unlist(OFF.list))
96}  # psv2magic
97
98
99