# R/common.R In smfishHmrf: Hidden Markov Random Field for Spatial Transcriptomic Data

#### Documented in checkErrorscheckErrorsMulticheckStopVerbosecheckStopVerboseMultifindDampFactorgetDeninitialIndicesinitialIndicesMultirelerrrelerr.multi

```#' @title checkErrors
#' @name checkErrors
#' @description check errors in mu and sigma (single variate)
#' @keywords internal
#check errors of input of tissue classification
checkErrors <- function(mu, sigma, err=NULL){
k <- length(mu)
if (k != length(sigma))
stop("The dimensions of 'mu' and 'sigma' do NOT match.")
if(! is.null(err))
if (any(err < 0)) stop("All 'err's have to be positive.")
}

#' @title getDen
#' @name getDen
#' @description calculate density given multivariate normal
#' @importFrom stats dnorm
#' @keywords internal
#Compute the density of all vertices.
getDen <- function(yunique, n.yunique, ymatch, mu, sigma){
k <- length(mu)
dy <- rep(yunique, each=k)
dmu <- rep(mu, n.yunique)
dsigma <- rep(sigma, n.yunique)
den <-  stats::dnorm(dy, mean = dmu, sd = dsigma)
den <- matrix(den, ncol=k, byrow=T)
den <- den[ymatch,]
den
}

#' @title relerr
#' @name relerr
#' @description compute the relative errors (single variate)
#' @keywords internal
#compute the relative errors
relerr <- function(x, y){
max(abs(x - y)) / (1 + max(abs(x), abs(y)))
}

#' @title initialIndices
#' @name initialIndices
#' @description initialize indices in the table (single variate)
#' @keywords internal
#Initialize indices
initialIndices <- function(y, nvert, mu, sigma, k){
index <- max.col(matrix(stats::dnorm(rep(y,k), mean=rep(mu, each=nvert),
sd=rep(sigma, each=nvert)), ncol=k))
indices <- do.call(cbind, lapply(1:k, function(x) indices==x))
rbind(indices, rep(0L,k))
}

#' @title checkStopVerbose
#' @name checkStopVerbose
#' @description check whether or not to stop (i.e. convergence)
#' @keywords internal
#check whether to stop the interation or not and whether ouput the curren state
checkStopVerbose <- function(muold, mu, sigmaold, sigma, err, it, maxit, verbose){
flag <- 0
remu <- relerr(muold, mu)
resigma <- relerr(sigmaold, sigma)
if (remu < err[1] && resigma < err[2] || it > maxit)
flag <- 1
if (verbose)
cat(paste("Iteration ", it, ": relative error of mean = ", signif(remu, 1),
"; sigma = ", signif(resigma, 1), "\n", sep=""))
flag
}

#' @title initialIndicesMulti
#' @name initialIndicesMulti
#' @description initialize indices for multivariate scenario
#' @importFrom pracma eps std
#' @keywords internal
#Initialize indices (multivariate), assume sub=FALSE, type="pure", subvox=NULL
initialIndicesMulti <- function(y, mu, sigma, k, dampFactor=NULL, forceDetectDamp=FALSE){
nvert <- dim(y)[1]
numdim <- dim(y)[2]

eps <- pracma::eps()
very_small <- mean(pracma::std(y)) * (1/k) * 0.0001

den = matrix(0, nrow=nvert, ncol=k)

for(j in 1:k){
determinant <- det(sigma[,,j])
invsigma <- array(0, dim(sigma[,,j]))
#choice 1=============
if(forceDetectDamp==FALSE & !is.null(dampFactor)){
sigma[,,j]<-sigma[,,j] + diag(numdim) * dampFactor[j]
invsigma <- solve(sigma[,,j])
determinant <- det(sigma[,,j])

}else if(forceDetectDamp==TRUE){
damp<-findDampFactor(sigma[,,j], factor=1.05, d_cutoff=1e-60, startValue=0.0001)
if(!is.null(damp)){
sigma[,,j]<-sigma[,,j] + diag(numdim) * damp
invsigma <- solve(sigma[,,j])
determinant <- det(sigma[,,j])
}
}
dist <- y - matrix(rep(t(mu[,j]),nvert),nrow=nvert,ncol=numdim,byrow=T)
exponent <- -0.5 * (rowSums((dist %*% invsigma) * dist))
const <- (numdim/2)*log(2*pi) + (1/2)*log(determinant)
exp_diff <- exponent - const
lim_exp_diff <- sapply(exp_diff, function(x) min(700, max(x, -700)))
den[,j] <- exp(lim_exp_diff)
}

indices <- max.col(den)

indices <- do.call(cbind, lapply(1:k, function(x) indices==x))
rbind(indices, rep(0L,k))
}

#' @title relerr.multi
#' @name relerr.multi
#' @description relative errors for multivariate scenario
#' @keywords internal
relerr.multi <- function(x, y){
max(abs(x - y)) / (1 + max(abs(x), abs(y)))
}

#' @title findDampFactor
#' @name findDampFactor
#' @description find dampening factor
#' @param sigma covariance matrix
#' @param factor step factor
#' @param d_cutoff determinant cutoff
#' @param startValue starting value to initialize the finding
#' @export
#' @examples
#' data(seqfishplus)
#' k<-dim(seqfishplus\$mu)[2]
#' damp<-array(0, c(k))
#' for(i in 1:k){
#'     di<-findDampFactor(seqfishplus\$sigma[,,i], factor=1.05, d_cutoff=1e-5, startValue=0.0001)
#'     damp[i]<-ifelse(is.null(di), 0, di)
#' }
findDampFactor <- function(sigma, factor=1.05, d_cutoff=1e-60, startValue=0.0001){
determinant <- det(sigma)
#factor <- 1.05
#d_cutoff<- 1e-60
numdim<-dim(sigma)[1]
damp<-NULL
if(determinant<d_cutoff){
print("Warning, need to add dampening factor to diagonal entries")
damp <- startValue
sigma_2 <- array(0, dim(sigma))
sigma_2 <- sigma + diag(numdim) * damp
determinant <- det(sigma_2)
if(determinant<d_cutoff){
repeat{
damp <- damp * factor
sigma_2 <- array(0, dim(sigma))
sigma_2 <- sigma + diag(numdim) * damp
determinant <- det(sigma_2)
if(determinant>d_cutoff) break
}
}else{
repeat{
damp <- damp / factor
sigma_2 <- array(0, dim(sigma))
sigma_2 <- sigma + diag(numdim) * damp
determinant <- det(sigma_2)
if(determinant<d_cutoff){
damp <- damp * factor
break
}
}
}
print(damp)
}
damp
}

#' @title checkStopVerboseMulti
#' @name checkStopVerboseMulti
#' @description check whether or not to stop (i.e. convergence) (multivariate)
#' @keywords internal
#check whether to stop the interation or not and whether ouput the curren state
checkStopVerboseMulti <- function(muold, mu, sigmaold, sigma, err, it, maxit, verbose){
flag <- 0
remu <- relerr.multi(muold, mu)
resigma <- relerr.multi(sigmaold, sigma)
if (remu < err[1] && resigma < err[2] || it > maxit)
flag <- 1
if (verbose)
cat(paste("Iteration ", it, ": relative error of mean = ", signif(remu, 1),
"; covariance = ", signif(resigma, 1), "\n", sep=""))
flag
}

#' @title checkErrorsMulti
#' @name checkErrorsMulti
#' @description check errors in mu and sigma (multivariate)
#' @keywords internal
#check errors of input of tissue classi
checkErrorsMulti <- function(mu, sigma, err=NULL){
k <- dim(mu)[2]
if (k != dim(sigma)[3])
stop("The dimensions of 'mu' and 'sigma' do NOT match.")
#if (!all(sigma >= 0))
#	stop("All 'sigma's have to be positive")
if(! is.null(err))
if (any(err < 0)) stop("All 'err's have to be positive.")
}
```

## Try the smfishHmrf package in your browser

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

smfishHmrf documentation built on Jan. 13, 2021, 12:54 p.m.