Nothing
rblockmult <- function(x,b){ # 2014-05-24 generalise to non-square b
m <- nrow(b)
n <- ncol(b)
stopifnot(ncol(x) %% m == 0) # 2014-05-24
r <- ncol(x)/m
res <- matrix(0, nrow = nrow(x), ncol = r * ncol(b))
x <- as.matrix(x) # in case they are from class 'Matrix' or other class
b <- as.matrix(b) # note that 'x' and 'b' cannot be vector, see nrow() and ncol() above.
for(i in 0:(r-1))
## 2022-10-19 r.h.s was: x[ , i*m + 1:m] %*% b
## but in Matrix v1.5-2 there is a change. This is from email by
## Mikael Jagan:
##
## "The reason is that the product of a pMatrix and a traditional matrix
## now a (row- or column-permuted) _M_atrix, rather than an equivalent
## _m_atrix, for consistency with other products involving _M_atrix."
##
## If 'b' is a permutation matrix (from package Matrix), then before the change
## the product on the r.h.s. below was an ordinary matrix (ok for the assignment)
## but after that it is a Matrix object and the assignment fails.
##
## This actually reveals a deficiency (or bug) in 'rblockmult'
## since it is clearly supposed to work with any matrix objects.
##
## was: res[ , i*n + 1:n] <- x[ , i*m + 1:m] %*% b
res[ , i*n + 1:n] <- x[ , i*m + 1:m] %*% b
res
}
## new: 2015-03-25
permute_var <- function(mat, perm = nrow(mat):1){ ## TODO: seems unused!
res <- mat[perm, perm]
if(mode(mat) == "numeric"){ # the above works for non-numeric matrices, as well.
P <- as(perm, "pMatrix") # lazy; todo: see also invPerm()
res2 <- P %*% mat %*% t(P)
stopifnot(all(res == res2))
}
res
}
permute_synch <- function(param, perm){
if(missing(perm)){
n <- .rnrow(param)
perm = n:1
}
P <- as(perm, "pMatrix") # lazy; todo: see also invPerm()
fu <- function(x){
if(is.list(x)){
for(i in seq(along = x)) # 2015-12-30 was: 1:length(x)
x[[i]] <- fu(x[[i]])
}else if(is.matrix(x))
x <- P %*% x %*% t(P)
else
x <- P %*% x # in case some components are vectors
x
}
fu(param)
}
.rnrow <- function(x){
if(is.list(x))
Recall(x[[1]])
else
NROW(x) # 2015-12-30 was: nrow(x)
}
## new: 2015-03-25
.ldl <- function(x){
R <- chol(x)
sqrtD <- diag(R)
d <- sqrtD^2 # lowercase d to emphasise that this is only the diagonal.
# todo: this may fail if x (or L) is (nearly) singular
R <- R / sqrtD # uses the recycling rule, i-th row of R is divided by sqrtD_i equivalent
# to dividing each row of L by sqrtD_i but slightly more convenient.
L <- t(R) # t() since chol() returns L'
diag(L) <- 1 # just to make sure that diagonal contains exact one's.
# theoretically (but not numerically), we have:
# stopifnot(all(x == L %*% diag(d) %*% t(L) ))
list(L = L, d = d)
}
.udu <- function(Sigma){
perm <- seq(nrow(Sigma), 1, by = -1) # 2015-12-30 was nros(Sigma):1, guard agains 0 rows
# P and t(P) are the same here, but for clarity use t(P) below
P <- as(perm, "pMatrix") # lazy
S <- t(P) %*% Sigma %*% P
wrk <- .ldl(S)
U <- P %*% wrk$L %*% t(P)
d <- wrk$d[perm]
## 2022-10-19 TODO: The note below now becomes urgent since a change in Matrix v1.5-2
## causes the result of matrix multiplication with permutation Matrix and 'matrix' be
## 'Matrix', not 'matrix' as before. This makes U of class Matrix.
##
## For now just converting U to matrix but redo this function as remarked below.
U <- as.matrix(U)
# TODO: the above is very lazy, could be done by simply permuting rows
# and columns. A simple check:
# D <- P %*% diag(wrk$d) %*% t(P)
# stopifnot(all(d == diag(D)))
list(U = U, d = d)
}
null_complement <- function(m, universe = NULL, na.allow = TRUE){
## edited 2015-07-10 to give error when both 'm' and 'universe' are NULL
if(na.allow && anyNA(m)){ # todo: this is probaly sensible only if all elem. of m are
# NA; could be refined, at least for the case when some
# columns of m are free from NA's
if(isNA(m)){
if(is.null(universe))
## Cannot determine the dimension of the space, so error.
stop("One of 'm' and 'universe' must be non-NULL.")
else
return(universe)
}##else m is assumed a matrix
if(is.null(universe))
universe <- diag(nrow = NROW(m))
if(all(is.na(m)))
res <- matrix(NA_real_, nrow = NROW(m), ncol = ncol(universe) - NCOL(m))
else{
## Zasega ostavyam kakto gornoto, vzh. komentara po-dolu.
##
res <- matrix(NA_real_, nrow = NROW(m), ncol = ncol(universe) - NCOL(m))
##
## TODO: rezultatat e lineyni komb. na kolonite na u2, ako vsyaka ot kolonite na
## 'm' e ili iztsyalo NA ili bez NA's. Inache (ako ima koloni s chisla i NA)
## tryabva oste rabota. PRI VSYAKO POLOZHENIE mi tryabva klas za parametrizirani
## pod-prostranstva, napr. e edin element za parametrite i vtori za bazisa.
## flags <- apply(m, 2, function(x) any(is.na(x)))
## wrk <- m[ , !flags] # select columns without any NA's
## u2 <- null_complement(wrk, universe = universe, na.allow = FALSE)
}
return(res)
}
if(is.null(universe))
return( Null(m) )
Null( cbind(m, Null(universe)) ) # compl of A w.r.t. B, where A is subspace of B
# equals complement of (A union B_orth)
} # for additional comments on orthogonal spaces see the comments at the end of this file.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.