R/types_util.R

In thirdwing/nimble: Flexible BUGS-compatible system for hierarchical statistical modeling and algorithm development

## adds explicit indexing to variables (using symtab), then expands indexing present on any variables
nl_expandNodeNames <- function(nodeNames, symtab, env) {
    nodeNames <- nl_addIndicesToVariables(nodeNames, symtab)
    nodeNames <- unlist(lapply(nodeNames, function(node) if(is.indexed(node)) nl_expandNodeIndex(node, env) else node))
    if(is.null(nodeNames))   return(character(0))
    return(nodeNames)
}


## Expands variables into their fully indexed form, e.g., 'y' is expanded to 'y[1:10]', using information in the symbolTable
nl_addIndicesToVariables <- function(nodeNames, symtab) {
    for(i in seq_along(nodeNames)) {
        nodeName <- nodeNames[i]
        varName <- nl_getVarNameFromNodeName(nodeName)
        if(!(varName %in% symtab$getSymbolNames()))   stop('variable not in symbol table')
        if(!is.indexed(nodeName) && (symtab$getSymbolField(varName, 'nDim') > 0)) {    ## nodeName has no indexing, and has dimension > 0
            maxs <- symtab$getSymbolField(varName, 'size')
            mins <- rep(1, length(maxs))
            indexStuff <- paste(mins, maxs, sep=':', collapse = ', ')
            nodeNames[i] <- paste0(varName, '[', indexStuff, ']')
        }
    }
    return(nodeNames)
}


## This is the same as nl_ExpandNodeIndex, except it takes a nodeExpr instead of a node char string
nl_expandNodeIndexExpr <- function(nodeExpr, env = parent.frame()) {
    if(length(nodeExpr)==1)  if(is.name(nodeExpr)) return(as.character(nodeExpr)) else stop('node expression with only one element, but not a variable name')
    indexExprs <- nodeExpr[-c(1,2)]
    indexStrs <- lapply(indexExprs, function(ind) as.character(eval(ind, envir=env)))
    numInd <- length(indexStrs)
    indexStrsCombined <- indexStrs[[numInd]]
    if(numInd > 1) for(i in seq(numInd-1, 1, by = -1)) {
        indexStrsCombined <- outer(indexStrs[[i]], indexStrsCombined, paste, sep=', ') }
    nodeStr <- as.character(nodeExpr[[2]])
    expandedNodes <- paste0(nodeStr, paste0('[',indexStrsCombined,']'))
    return(expandedNodes)
}

## same as nl_vectorizedExpandNodeIndex but takes nodeExprs instead of nodes
nl_vectorizedExpandNodeIndexExprs <- function(nodeExprs, env = parent.frame()) {
    return(unlist(lapply(nodeExprs, function(n) nl_expandNodeIndexExpr(n, env))))
}

## Expands the indexing of a single node name string, e.g., 'x[1:3]' is expanded to c('x[1]', 'x[2]', 'x[3]')
nl_expandNodeIndex <- function(node, env = parent.frame()) {
    nodeExpr <- parse(text=node, keep.source = FALSE)[[1]]
    if(length(nodeExpr)==1)  if(is.name(nodeExpr)) return(as.character(nodeExpr)) else stop('node expression with only one element, but not a variable name')
    indexExprs <- nodeExpr[-c(1,2)]
    indexStrs <- lapply(indexExprs, function(ind) as.character(eval(ind, envir=env)))
    numInd <- length(indexStrs)
    indexStrsCombined <- indexStrs[[numInd]]
    if(numInd > 1) for(i in seq(numInd-1, 1, by = -1)) {
        indexStrsCombined <- outer(indexStrs[[i]], indexStrsCombined, paste, sep=', ') }
    nodeStr <- as.character(nodeExpr[[2]])
    expandedNodes <- paste0(nodeStr, paste0('[',indexStrsCombined,']'))
    return(expandedNodes)
}


nl_vectorizedExpandNodeIndex <- function(nodes, env = parent.frame()) {
    return(unlist(lapply(nodes, function(n) nl_expandNodeIndex(n, env))))
}


## checks that all nodeNames are present in model
nl_checkNodeNamesInModel <- function(model, nodeNames, determOnly = FALSE, stochOnly = FALSE) {
    if(!all(nodeNames %in% model$getNodeNames()))      stop('all node names not in model')
    if(determOnly) if(!all(nodeNames %in% model$getNodeNames(determOnly = TRUE)))      stop('all node names are not deterministic')
    if(stochOnly)  if(!all(nodeNames %in% model$getNodeNames(stochOnly = TRUE)))       stop('all node names are not stochastic')
}

## checks that all varNames are present in model
nl_checkVarNamesInModel <- function(model, varNames) {
    if(!all(varNames %in% model$getVarNames(logProb = TRUE)))      stop('all variable names are not in model')
}

#nl_nodeVectorReadyNodes <- function(model, nodeNames, includeData = TRUE){
#	expandedNames <- model$expandNodeNames(nodeNames)
#	sortedNames = model$topologicallySortNodes(expandedNames)
#	if(!includeData){
#		keepNodes = rep(NA, length(sortedNames) ) 
#		for(i in seq_along(keepNodes) ) 
#			keepNodes[i] <- !(model$isData(sortedNames[i]) )
#		sortedNames <- sortedNames[keepNodes]
#	}
#	return(sortedNames)
#}

nl_nodeVectorReadyNodes <- function(model, nodeNames, includeData = TRUE){
	GIDs <- model$modelDef$nodeName2GraphIDs(nodeNames)
	sortedIDs <- sort(GIDs)
	if(includeData == FALSE)
		sortedIDs = sortedIDs[!model$isDataFromGraphID(GIDs)]
	return(model$modelDef$maps$graphID_2_nodeFunctionName[sortedIDs])
}




## returns a list of variable names, and flat index ranges
nl_createVarsAndFlatIndexRanges <- function(nodeNames, symtab) {
    varAndIndexes <- lapply(nodeNames, function(nn) list(var = nl_getVarNameFromNodeName(nn), ind = if(is.indexed(nn)) unlist(as.list(parse(text=nn, keep.source = FALSE)[[1]][-c(1,2)])) else numeric(0)))
    varAndFlatIndexes <- lapply(varAndIndexes, function(vai) list(var = vai$var, flatIndex = nl_determineFlatIndex(ind=vai$ind, maxs=symtab$getSymbolField(vai$var, 'size'))))
    listsOfIndexes <- list()
    for(vafi in varAndFlatIndexes) listsOfIndexes[[vafi$var]] <- list(var = vafi$var, inds = c(listsOfIndexes[[vafi$var]]$inds, vafi$flatIndex))
    listsOfAggIndexes <- lapply(listsOfIndexes, function(loi) list(var=loi$var, agInd=nl_aggregateConsecutiveBlocks(loi$inds)))
    separatedListsOfAggIndexes <- list()
    for(loai in listsOfAggIndexes) for(ind in loai$agInd) separatedListsOfAggIndexes <- c(separatedListsOfAggIndexes, list(list(var=loai$var, ind=ind)))
    return(separatedListsOfAggIndexes)
}


## determines the (scalar) flatIndex, from numeric vectors of indicies, and max_indicies
nl_determineFlatIndex = function(ind, maxs) {
		
    if(length(ind) != length(maxs))   stop('number of indices and dimensions are different')
    if(any(ind > maxs))               stop('some indices exceed dimensions')
    if(any(ind < 1))                  stop('non-positive index')
    nDim <- length(ind)
    if(nDim==0) return(1)
    if(nDim==1) return(ind)
    max_prods <- 1
    for(iDim in 2:nDim) max_prods[iDim] <- prod(maxs[1:(iDim-1)])
    indShifted <- c(ind[1], ind[-1]-1)
    flatIndex <- sum(indShifted * max_prods)
    return(flatIndex)
}


## takes a vector of flat indices, e.g. c(1,2,3,6,7,100),
## returns a list of pairs: (index_start, index_finish), e.g. [[1]] c(1,3)  [[2]] c(6,7)  [[3]] c(100,100)
## If anyone has an implementation which beats O(n) linear time, then let's use it.
nl_aggregateConsecutiveBlocks <- function(ind) {
    if(length(ind) == 0)     return(list())
#    ind <- unique(sort(ind))
    indDiffLogical <- c(TRUE, diff(ind) != 1)
    aggregated <- list()
    for(i in seq_along(indDiffLogical)) {
        if(indDiffLogical[i])     aggregated[[length(aggregated)+1]] <- c(ind[i], ind[i])
        if(!indDiffLogical[i])    aggregated[[length(aggregated)]][2] <- ind[i]
    }
    return(aggregated)
}

nl_removeNodeNamesNotInSymbolTable <- function(nodeNames, st) {
    return(nodeNames[nl_getVarNameFromNodeName(nodeNames) %in% st$getSymbolNames()])
}


nl_getVarNameFromNodeName <- function(nodeName)    gsub('\\[.*', '', nodeName)

#nimDim <- function(object){
#	rDim = dim(object)
#	if(is.null(rDim) ) {
#		if(is.vector(object) ) 
#			return(length(object) ) 
#		stop('dim called on object ', as.character(substitute(object) ), ' for which dim is undefined')
#	}
#	return(dim(object) ) 
#}

expandMVNames <- function(mv, varNames){
	sizeList = mv$sizes
#	varNames = names(sizeList)
	nodeNames = NA
	nodeIndex = 0 
	for(i in seq_along(varNames) ){
		baseName = varNames[i]
		dims = mv$symTab$getSymbolObject(baseName)$size
		if(length(dims) == 0){
			nodeNames[nodeIndex+1] = baseName
			nodeIndex <- nodeIndex + 1
			}
		else{
            mins <- rep(1, length(dims))
            indexStuff <- paste(mins, dims, sep=':', collapse = ', ')
            compactNodeNames <- paste0(baseName, '[', indexStuff, ']')
            expandedNodeNames <- nl_expandNodeIndex(compactNodeNames)
            nodeNames[nodeIndex + 1:length(expandedNodeNames)] <- expandedNodeNames
            nodeIndex = nodeIndex + length(expandedNodeNames)
		}
	}
	return(nodeNames)
}

as.matrix.modelValuesBaseClass <- function(mv, varNames){
	if(missing(varNames))
			varNames <- mv$varNames
	nrows = getsize(mv)
	flatNames = expandMVNames(mv, varNames)
	ans <- matrix(0.1, nrow = nrows, ncol = length(flatNames))
	colIndex = 1
	for(i in seq_along(varNames)){
		#ans[1:nrows, colIndex + 1:prod(mv$sizes[[varNames[i] ]]) ] <- modelValuesElement2Matrix(mv, varNames[i])
		.Call('fastMatrixInsert', ans, modelValuesElement2Matrix(mv, varNames[i]) , as.integer(1), as.integer(colIndex) ) 		
		#This method is faster than the commented out method, especially if there are a lot of variables (as apposed to nodes)
		colIndex = colIndex + prod(mv$sizes[[varNames[i] ]])
		}
	colnames(ans) <- flatNames
	return(ans)

}

as.matrix.CmodelValues <- function(mv, varNames){
	if(missing(varNames))
			varNames <- mv$varNames
	nrows = getsize(mv)
	flatNames = expandMVNames(mv, varNames)
	ans <- matrix(0.1, nrow = nrows, ncol = length(flatNames))
	colIndex = 1
	for(i in seq_along(varNames)){
		#ans[1:nrows, colIndex + 1:prod(mv$sizes[[varNames[i] ]]) ] <- modelValuesElement2Matrix(mv, varNames[i])
		.Call('fastMatrixInsert', ans, modelValuesElement2Matrix(mv, varNames[i]) , as.integer(1), as.integer(colIndex) ) 		
		#This method is faster than the commented out method, especially if there are a lot of variables (as apposed to nodes)
		colIndex = colIndex + prod(mv$sizes[[varNames[i] ]])
		}
	colnames(ans) <- flatNames
	return(ans)
}


modelValuesElement2Matrix <- function(mv, varName){
	if(length(varName) != 1)
		stop('modelValuesElement2Matrix is a call for a single variable. For multiple variables, use modelValues2Matrix')
	matrix( as.numeric(unlist(mv[[varName]]) ), ncol = prod(mv$sizes[[varName]]), byrow = TRUE )
}


Cmatrix2mvOneVar <- function(mat, mv, varName, k){
	ptr <- mv$componentExtptrs[[varName]]
	if(inherits(ptr, 'externalptr'))
		.Call('matrix2VecNimArr', ptr, mat, rowStart = as.integer(1), rowEnd = k )
	else
		stop('varName not found in modelValues')
}

Rmatrix2mvOneVar <- function(mat, mv, varName, k){
	if( mv$symTab$symbols[[varName]][['type']] == 'double'){
		storage.mode(mat) <- 'double'
		len <- ncol(mat)
		.Call('matrix2ListDouble', mat, mv[[varName]], listStartIndex = as.integer(1), RnRows = k, Rlength = as.integer(mv$sizes[[varName]]) )
	}
	if( mv$symTab$symbols[[varName]][['type']] == 'int'){
		storage.mode(mat) <- 'integer'
		len <- ncol(mat)
		.Call('matrix2ListInt', mat, mv[[varName]], listStartIndex = as.integer(1), RnRows = k, Rlength = as.integer(mv$sizes[[varName]]) )
	}
}

matrix2mv <- function(mat, mv){
	k <- nrow(mat)
	if(mv$getSize() < k)
		mv$resize(k)
	mvVarNames <- mv$varNames
	mvSizes <- mv$sizes
	colNames <- colnames(mat)
	colNames <- removeIndexing(colNames)
	uniqueColNames <- unique(colNames)
	if(inherits(mv, 'CmodelValues')	){
		for(vN in uniqueColNames){
			totVals <- prod(mvSizes[[vN]])
			varInds <- colNames == vN
			if(totVals != sum(varInds) )
				stop('matrix2mv halted because dimensions of variables do not match')
			subMatrix <- as.matrix(mat[, varInds])
			Cmatrix2mvOneVar(subMatrix, mv, vN, k)
		}
	}
	else if(inherits(mv, 'modelValuesBaseClass') ){
		for(vN in uniqueColNames){
			totVals <- prod(mvSizes[[vN]])
			varInds <- colNames == vN
			if(totVals != sum(varInds) )
				stop('matrix2mv halted because dimensions of variables do not match')
			subMatrix <- as.matrix(mat[, varInds])
			Rmatrix2mvOneVar(subMatrix, mv, vN, k)
		}
	}
	else
		stop('argument mv is neither a CmodelValues or RmodelValues object')
}

#modelValues2Matrix <-function(mv, varNames){
#	if(missing(varNames))
#			varNames <- mv$varNames
#	nrows = getsize(mv)
#	flatNames = expandMVNames(mv, varNames)
#	ans <- matrix(0.1, nrow = nrows, ncol = length(flatNames))
#	colIndex = 1
#	for(i in seq_along(varNames)){
#		#ans[1:nrows, colIndex + 1:prod(mv$sizes[[varNames[i] ]]) ] <- modelValuesElement2Matrix(mv, varNames[i])
#		.Call('fastMatrixInsert', ans, modelValuesElement2Matrix(mv, varNames[i]) , as.integer(1), as.integer(colIndex) ) 		
#		#This method is faster than the commented out method, especially if there are a lot of variables (as apposed to nodes)
#		colIndex = colIndex + prod(mv$sizes[[varNames[i] ]])
#		}
#	colnames(ans) <- flatNames
#	return(ans)
#}