# R/utilities_grips.R In gRim: Graphical Interaction Models

#### Documented in AIC.gips_fit_classBIC.gips_fit_classconcentrationconcentration.gips_fit_classfast_covgenerate_n01glance.gips_fit_classimpose_zerologLik.gips_fit_classprint.gips_fit_classsigma.gips_fit_classsummary.gips_fit_class

```#' Fast computation of covariance / correlation matrix
#'
#' @param x a numeric matrix(like object).
#' @param center,scale Should columns in x be centered and/or scaled
#' @export
fast_cov <- function(x, center=TRUE, scale=TRUE){
A <- scale(x, center=center, scale=scale) |> crossprod()
return(A / (nrow(x) - 1))
}

#' Genrate matrix of N(0, 1) variables
#'
#' @param n.obs Number of rows
#' @param nvar Number of columns
#' @param seed Seed for random number generator
#'
#' @export
generate_n01 <- function(n.obs, nvar, seed=2022){
set.seed(seed)
d <-matrix(rnorm(n.obs * nvar), nrow=n.obs)
d
}

solve_fun <- function(A, b=NULL){
solve_qr <- function(A, b=NULL){
if (is.null(b))
solve.qr(qr(A))
else
solve.qr(qr(A)) %*% b
}

## solve_ginv <- function(A, b=NULL){
##     if (is.null(b))
##         MASS::ginv(A)
##     else
##         MASS::ginv(A) %*% b
## }

## solveM <- function(A, b=NULL){
##     if (is.null(b))
##         Matrix::solve(A)
##     else
##         Matrix::solve(A, b)
## }

solve_qr(A, b)
}

ggm_logL <- function(S, K, nobs){
nobs * (log(det(K)) - sum(S * K)) / 2
}

## UTILITIES

.get.time <- function(){
Sys.time()
}

.get.diff.time <- function(t0, units="secs"){
unit_ <- match.arg(units, c("secs", "millisecs"))
d <- as.numeric(difftime(Sys.time(), t0, units="secs"))
if (identical(units, "millisecs"))
d <- d * 1000
d
}

matchKSigma <- function(object){
if (!inherits(object, "gips_fit_class"))
stop("'object' is not a gips_git_class object")
## sum(abs((object\$K %*% object\$Sigma) - diag(1, nrow(object\$K))))
a <- (object\$K %*% object\$Sigma) - diag(1, nrow(object\$K))
## sqrt(sum(a^2))
norm(a, type="F")
}

## #' Parameter equality
## #'
## #' @param p2a,p2b Parameter lists.
## #' @param eps A small number.
## #'
## #' @export
## #' @rdname utilities
## parmEq <- function(p2a, p2b, eps=1e-4){
##   max(abs(p2a\$K - p2b\$K)) / max(abs(p2a\$K)) < eps && max(abs(p2a\$Sigma - p2b\$Sigma)) / max(abs(p2a\$Sigma)) < eps
## }

#' Impose zeros in matrix entries which do not correspond to an edge.
#'
#' @param emat Edge matrix (2 x p matrix)
#' @param K Matrix; typically a concentration matrix.
#'
#' @export
impose_zero <- function(emat, K){
emc <- as_emat_complement(emat, nrow(K))
if (ncol(emc) == 0){
return(K)
}
emc <- t.default(emc)
K[rbind(emc, emc[2:1,])]  <- 0
return(K)
}

### Internal .functions

.sol <- function(A){ ## FIXME : silly name
.solve2x2 <- function(A){
##d <- A[1, 1] * A[2, 2] - A[1, 2] * A[2, 1] ## FIXME: Handle d=0 case
d <- A[1] * A[4] - A[2] * A[3] ## FIXME: Handle d=0 case
#B <- A
B <- c(0, 0, 0, 0)
B[1] <- A[4]
B[2] <- -A[2]
B[3] <- -A[3]
B[4] <- A[1]
B <- B * (1 / d)
dim(B) = c(2L, 2L)
B
}

if (dim(A)[1] == 2) .solve2x2(A)
else solve.default(A) ## FIXME: Could do chol2inv(chol(A))
}

.checkS <- function(S){
if (!is.matrix(S))
stop("Second argument is not a matrix!")

if (dim(S)[1] != dim(S)[2])
stop("Second argument is not a square matrix!")
}

.penalizeS <- function(S, lambda=0){
S + diag(lambda, nrow(S))
}

## Just once
.make_rc <- function(glist, S){
if (!is.matrix(glist))
glist <- do.call(cbind, glist)
vn    <- unique.default(sort(glist))
amat  <- matrix(0L, nrow=length(vn), ncol=length(vn))
amat[t.default(glist)] <- 1L
amat <- amat + t.default(amat)
dimnames(amat) <- dimnames(S)
lt   <- lower.tri(amat)
rrr  <- which.arr.index(lt * amat != 0)
colnames(rrr) <- c("row", "col")
rrr
}

#' @title Utilities for gRips
#'
#' @name utilities_grips

#' @export
#' @rdname utilities_grips
#' @param object Model object.
#' @param ... Additional arguments; currently not used.
#' @param k Penalty parameter for calculating AIC; only k=2 gives genuine AIC.
#' @param x Object to be printed.
logLik.gips_fit_class <- function(object, ...){

nobs <- nobs(object)# unname(object\$details["nobs"])
if (is.na(nobs))
stop("'nobs' not given; can not compute log L\n")

trKS  <- unname(object\$details\$trKS)
nparm <- unname(object\$details\$dim)

##str(list(nobs, trKS))
out  <- c(nobs * (log(det(object\$K)) - trKS) / 2)

## FIXME: Perhaps need number of parameters in specified model.
attr(out, "nobs")  = nobs
attr(out, "nparm") = nparm
attr(out, "df")    = nparm

class(out) <- "logLik"
out
}

#' @method AIC gips_fit_class
#' @export
#' @rdname utilities_grips
AIC.gips_fit_class <- function(object, ..., k=2){
ll <- logLik(object)
-2 * as.numeric(ll) + k * attr(ll, "nparm")
}

#' @method BIC gips_fit_class
#' @export
#' @rdname utilities_grips
BIC.gips_fit_class <- function(object, ...){
ll <- logLik(object)
-2 * as.numeric(ll) + log(attr(ll, "nobs")) * attr(ll, "nparm")
}

#' @method sigma gips_fit_class
#' @export
#' @rdname utilities_grips
sigma.gips_fit_class <- function(object, ...) {
object\$Sigma
}

#' @export
#' @rdname utilities_grips
concentration <- function(object, ...) {
UseMethod("concentration")
}

#' @method concentration gips_fit_class
#' @export
#' @rdname utilities_grips
concentration.gips_fit_class <- function(object, ...) {
object\$K
}

#' @export
#' @rdname utilities_grips
print.gips_fit_class <- function(x, ...){
## cat("Method: ", x\$method, "\n")
## xx <- c(method=x\$method, eng=x\$engine, x\$details)
## print(as.data.frame(x\$details))
xx <- x\$details
xx <- as.data.frame(xx)
xx\$engine <- NULL
xx\$ncore <- NULL
print(as.data.frame(xx))
xx
}

#' @export
#' @rdname utilities_grips
summary.gips_fit_class <- function(object, ...) {
xx <- as.data.frame(object\$details)
xx\$engine <- NULL
xx\$ncore <- NULL
## xx\$tpi <- with(xx, time / iter)
## xx\$tpe <- with(xx, time / (dim-idim))
xx
}

#' @export
#' @rdname utilities_grips
glance.gips_fit_class <- function(x, ...) {
as.data.frame(x\$details)
}
```

## Try the gRim package in your browser

Any scripts or data that you put into this service are public.

gRim documentation built on Sept. 11, 2024, 7:02 p.m.