Nothing
R/learnDbnStruct3dParDeg1.R
In TGS: Rapid Reconstruction of Time-Varying Gene Regulatory Networks
Defines functions
LearnDbnStructMo1Clr3Ser
learnDbnStructMo1Layer3dParDeg1_v2
learnDbnStructMo1Layer3dParDeg1
learnDbnStructLayer3dParDeg1
# #Goal: Infer Dynamic Bayesian Network (DBN)
#
# #---------------------------------
# # Begin: Step 2: Decomposing the network
# #---------------------------------
#
#
# #' Unrolled DBN structure learning with all possible Markov Orders.
# #'
# #' Unrolled DBN structure learning with all possible Markov Orders.
# #' Candidate parents: The target node itself and its CLR net neighbours at any previous and current time pt.
# #'
# #' @import bnstruct
# #' @import ggm
# #' @import foreach
# #' @import doParallel
# #'
# #' @param input.data.discr.3D 3D input data Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
# #' @param mi.net.adj.matrix Adjacency matrix of the mutual information network. Rownames and colnames should be node names.
# #' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized
# #' @param num.nodes number of nodes
# #' @param num.timepts number of timepoints
# #'
# #' @return Unrolled DBN adjacency matrix
# #'
# #' @export
# learnDbnStruct3dParDeg1 <- function(input.data.discr.3D, mi.net.adj.matrix, num.discr.levels, num.nodes, num.timepts)
# {
# # load('dream3.yeast1.size10.trajectory.3D.Rdata')
# # input.data.3D <- dream3.yeast1.size10.trajectory.3D.data
# # load('mi.net.adj.matrix.Rdata')
#
# # In 'di.net.adj.matrix', rows are src nodes and cols are tgt nodes
# di.net.adj.matrix <- mi.net.adj.matrix
# di.net.adj.matrix[1:base::nrow(di.net.adj.matrix), 1:base::ncol(di.net.adj.matrix)] <- 0
#
# no_cores <- num.nodes
# cl <- parallel::makeCluster(no_cores)
# doParallel::registerDoParallel(cl)
#
# #base::print('going inside foreach')
#
# # For each central node, learn local DBN struct in parallel.
# # Package names, specified by '.packages' must be copied to each worker core.
# # Current environmental variables that are referenced inside the foreach loop are copied
# # to each worker automatically.
# # local.unrolled.DBN.adj.matrix.list <- foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = 'bnstruct') %dopar%
# # local.unrolled.DBN.adj.matrix.list <- as.list(1:10)
#
# # local.unrolled.DBN.adj.matrix.list <- foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = 'bnstruct') %:% when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %dopar%
# # {
# # print('now inside foreach 1')
# #
# # centralNodeName <- rownames(mi.net.adj.matrix)[centralNodeIdx]
# #
# # # List names of the central node's neighbours in mi.net.adj.matrix
# # nbghNames <- c()
# # if (sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
# # {
# # for (nbrIdx in 1:n)
# # {
# # if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
# # {
# # nbghNames <- rownames(mi.net.adj.matrix)[nbrIdx]
# # break
# # }
# # }
# # }
# # else if (sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
# # {
# # nbghNames <- rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
# # }
# #
# # local.net.node.names <- c(centralNodeName, nbghNames)
# #
# # local.input.data.3D <- input.data.discr.3D[, local.net.node.names, ]
# #
# # #---------------------------------
# # # Begin: Local Unrolled DBN struct learning
# # #---------------------------------
# #
# # sampleIdx <- 1
# #
# # local.DBN.input.data <- matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
# # ncol = dim(local.input.data.3D)[2])
# #
# # for (time.pt in 2:dim(local.input.data.3D)[1])
# # {
# # data.to.combine <- matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# # ncol = dim(local.input.data.3D)[2])
# # local.DBN.input.data <- cbind(local.DBN.input.data, data.to.combine)
# # }
# #
# # local.DBN.input.data.var.names <- c()
# # for (time.pt in 1:dim(local.input.data.3D)[1]) {
# # for (var.name in dimnames(local.input.data.3D)[2]) {
# # local.DBN.input.data.var.names <- c(local.DBN.input.data.var.names,
# # paste(var.name, as.character(time.pt), sep = "_t"))
# # }
# # }
# # colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
# #
# # if (dim(local.input.data.3D)[3] > 1) # If there are multiple samples per time pt
# # {
# # for (sampleIdx in 2:dim(local.input.data.3D)[3])
# # {
# # sample.to.combine <- matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
# # ncol = dim(local.input.data.3D)[2])
# #
# # for (time.pt in 2:dim(local.input.data.3D)[1])
# # {
# # data.to.combine <- matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# # ncol = dim(local.input.data.3D)[2])
# # sample.to.combine <- cbind(sample.to.combine, data.to.combine)
# # }
# #
# # local.DBN.input.data <- rbind(local.DBN.input.data, sample.to.combine)
# # }
# # }
# #
# # rownames(local.DBN.input.data) <- as.vector(unlist(dimnames(local.input.data.3D)[3]))
# #
# # local.unrolled.DBN.adj.matrix <- matrix(0, nrow = ncol(local.DBN.input.data), ncol = ncol(local.DBN.input.data),
# # dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# #
# # # 'fixed = FALSE' represents that the given pattern is a regular expression.
# # local.unrolled.DBN.central.node.indices <- grep(paste('^', centralNodeName, sep = ''),
# # colnames(local.unrolled.DBN.adj.matrix),
# # fixed = FALSE)
# #
# # # Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# # local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
# #
# # #---------------------------------
# # # End: Local Unrolled DBN struct learning
# # #---------------------------------
# #
# # local.unrolled.DBN.adj.submatrix
# # }
# #
# # print(local.unrolled.DBN.adj.matrix.list)
#
# local.unrolled.DBN.adj.matrix.list <- foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = 'bnstruct') %:% when(base::sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %dopar%
# {
#
# # print('now inside foreach')
#
# # centralNodeIdx <- 1
#
# centralNodeName <- base::rownames(mi.net.adj.matrix)[centralNodeIdx]
#
# # if central node does not have any neighbour in mutual information net
# if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 0)
# {
# next
# }
#
# # List names of the central node's neighbours in mi.net.adj.matrix
# nbghNames <- base::c()
# if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
# {
# for (nbrIdx in 1:n)
# {
# if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
# {
# nbghNames <- base::rownames(mi.net.adj.matrix)[nbrIdx]
# break
# }
# }
# }
# else if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
# {
# nbghNames <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
# }
#
# local.net.node.names <- base::c(centralNodeName, nbghNames)
#
# local.input.data.3D <- input.data.discr.3D[, local.net.node.names, ]
#
#
# # #---------------------------------
# # # Begin: Local BN struct learning
# # #---------------------------------
# # local.BN.input.data <- bnstruct::BNDataset(local.discr.data,
# # variables = colnames(local.discr.data),
# # discreteness = rep(TRUE, ncol(local.discr.data)),
# # node.sizes = rep(2, ncol(local.discr.data)),
# # starts.from = 0)
# #
# # # Default params: scoring.func = "BDeu"
# # local.BN <- bnstruct::learn.network(local.discr.data, algo = 'sm')
# # plot(local.BN)
# #
# # # Rows are src nodes and cols are tgt nodes
# # local.BN.adj.matrix <- bnstruct::dag(local.BN)
# # local.BN.adj.matrix <- matrix(local.BN.adj.matrix, nrow = length(local.BN@variables),
# # ncol = length(local.BN@variables),
# # dimnames = c(list(local.BN@variables), list(local.BN@variables)))
# #
# # local.BN.adj.matrix[, centralNodeName]
# # #---------------------------------
# # # End: Local BN struct learning
# # #---------------------------------
#
#
# #---------------------------------
# # Begin: Local Unrolled DBN struct learning
# #---------------------------------
#
# sampleIdx <- 1
#
# local.DBN.input.data <- base::matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
# ncol = base::dim(local.input.data.3D)[2])
#
# for (time.pt in 2:base::dim(local.input.data.3D)[1])
# {
# data.to.combine <- base::matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# ncol = base::dim(local.input.data.3D)[2])
# local.DBN.input.data <- base::cbind(local.DBN.input.data, data.to.combine)
# }
#
# local.DBN.input.data.var.names <- base::c()
# for (time.pt in 1:base::dim(local.input.data.3D)[1]) {
# for (var.name in base::dimnames(local.input.data.3D)[2]) {
# local.DBN.input.data.var.names <- base::c(local.DBN.input.data.var.names,
# base::paste(var.name, as.character(time.pt), sep = "_t"))
# }
# }
# base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
#
# if (base::dim(local.input.data.3D)[3] > 1) # If there are multiple samples per time pt
# {
# for (sampleIdx in 2:base::dim(local.input.data.3D)[3])
# {
# sample.to.combine <- base::matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
# ncol = base::dim(local.input.data.3D)[2])
#
# for (time.pt in 2:base::dim(local.input.data.3D)[1])
# {
# data.to.combine <- base::matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# ncol = base::dim(local.input.data.3D)[2])
# sample.to.combine <- base::cbind(sample.to.combine, data.to.combine)
# }
#
# # local.DBN.input.data.var.names <- c()
# # for (time.pt in 1:dim(local.input.data.3D)[1]) {
# # for (var.name in dimnames(local.input.data.3D)[2]) {
# # local.DBN.input.data.var.names <- c(local.DBN.input.data.var.names,
# # paste(var.name, as.character(time.pt), sep = "_t"))
# # }
# # }
# # colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
#
# local.DBN.input.data <- base::rbind(local.DBN.input.data, sample.to.combine)
# }
# }
#
# base::rownames(local.DBN.input.data) <- base::as.vector(base::unlist(base::dimnames(local.input.data.3D)[3]))
#
# # local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
# # discreteness = rep(TRUE, ncol(local.DBN.input.data)),
# # variables = colnames(local.DBN.input.data),
# # node.sizes = rep(num.discr.levels, ncol(local.DBN.input.data)),
# # starts.from = 0,
# # num.time.steps = nrow(local.input.data.3D))
#
#
# # Begin: Uncomment this section after testing
# local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
# discreteness = base::rep(TRUE, base::ncol(local.DBN.input.data)),
# variables = base::colnames(local.DBN.input.data),
# node.sizes = base::rep(num.discr.levels, base::ncol(local.DBN.input.data)),
# num.time.steps = base::nrow(local.input.data.3D))
#
# # algo = "mmhc", scoring.func = "BDeu", no layering
# # local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# # num.time.steps =
# # bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
#
# # algo = "mmhc", scoring.func = "BDeu", with layering.
# # A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- base::c()
# for (time.pt in 1:num.timepts)
# {
# layers <- base::c(layers, base::rep(time.pt, base::dim(local.input.data.3D)[2]))
# }
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# initial.network = 'random.chain', seed = 12345,
# algo = 'mmhc',
# scoring.func = 'BDeu',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
#
#
# # algo = "sm", scoring.func = "BDeu", no layering
# # local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# # algo = 'sm',
# # num.time.steps =
# # bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
#
#
# # # The following four lines must be executed at the same time
# # save.plot.to.filename = paste(paste('LocalUnrolledDbn', centralNodeName, sep = '_'), '.jpg', sep = '')
# # jpeg(file = paste('LocalUnrolledDbn_', centralNodeName, '.jpg', sep = ''))
# # plot(local.unrolled.DBN)
# # dev.off()
#
# # Extracting the adjacency matrix of the local DBN
# local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
# local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
# nrow = base::length(local.unrolled.DBN@variables),
# ncol = base::length(local.unrolled.DBN@variables),
# dimnames = base::c(base::list(local.unrolled.DBN@variables), base::list(local.unrolled.DBN@variables)))
# # End: Uncomment this section after testing
#
# # # Begin: This section is for testing
# # local.unrolled.DBN.adj.matrix <- matrix(0,
# # nrow = ncol(local.DBN.input.data),
# # ncol = ncol(local.DBN.input.data),
# # dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# # # End: This section is for testing
#
# # 'fixed = FALSE' represents that the given pattern is a regular expression.
# local.unrolled.DBN.central.node.indices <- base::grep(base::paste('^', centralNodeName, sep = ''),
# base::colnames(local.unrolled.DBN.adj.matrix),
# fixed = FALSE)
#
# # Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
#
# #---------------------------------
# # End: Local Unrolled DBN struct learning
# #---------------------------------
#
# # Return value for each 'foreach' iteration
# local.unrolled.DBN.adj.submatrix
# }
#
# # print('foreach ended')
#
# # Shut down the cluster
# parallel::stopCluster(cl)
#
# # print('after stopCluster')
#
# # Begin: Unrolled DBN struct learning
#
# # Begin: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
# unrolled.DBN.adj.matrix.node.names <- base::c()
# for (time.pt in 1:num.timepts) {
# for (node.name in base::dimnames(input.data.discr.3D)[2]) {
# unrolled.DBN.adj.matrix.node.names <- base::c(unrolled.DBN.adj.matrix.node.names,
# base::paste(node.name, as.character(time.pt), sep = "_t"))
# }
# }
#
# unrolled.DBN.adj.matrix <- base::matrix(0, nrow = base::length(unrolled.DBN.adj.matrix.node.names),
# ncol = base::length(unrolled.DBN.adj.matrix.node.names),
# dimnames = base::c(base::list(unrolled.DBN.adj.matrix.node.names),
# base::list(unrolled.DBN.adj.matrix.node.names)))
#
# # End: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
#
# for (list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list))
# {
# # 'submatrix.to.combine'
# submatrix.to.combine <- local.unrolled.DBN.adj.matrix.list[[list.idx]]
#
# unrolled.DBN.adj.matrix[base::rownames(submatrix.to.combine), base::colnames(submatrix.to.combine)] <- submatrix.to.combine
# }
#
# # End: Unrolled DBN struct learning
#
# # Return the final directed net of the current iteration
# return (unrolled.DBN.adj.matrix)
# }
#' Unrolled DBN structure learning with Markov Order 0 and 1.
#'
#' Unrolled DBN structure learning with Markov Order 0 and 1.
#' Candidate parents: The target node itself and its CLR net neighbours at immediately previous and current time pt.
#'
#' @param input.data.discr.3D Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
#' @param mi.net.adj.matrix Adjacency matrix of the mutual information network. Rownames and colnames should be node names.
#' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized for performing the DBN structure learning.
#' @param num.nodes number of nodes
#' @param num.timepts number of timepoints
#' @param output.dirname output directory to store files
#'
#' @return Unrolled DBN adjacency matrix
#'
#' @keywords internal
#' @noRd
#' @importFrom foreach %do% %:% when
learnDbnStructLayer3dParDeg1 <- function(input.data.discr.3D, mi.net.adj.matrix, num.discr.levels, num.nodes, num.timepts, output.dirname="./OUTPUT") {
if(!base::is.array(input.data.discr.3D))
{
base::stop("Error in learnDbnStructLayer3dParDeg1 input.data.discr.3D is not an array")
}
if(!base::is.matrix(mi.net.adj.matrix))
{
base::stop("Error in learnDbnStructLayer3dParDeg1 mi.net.adj.matrix is not a matrix")
}
num.time.trans <- (num.timepts - 1)
no_cores <- num.nodes
cl <- parallel::makeCluster(no_cores, outfile = base::paste(output.dirname, 'outfile.txt', sep = '/' ))
doParallel::registerDoParallel(cl)
centralNodeIdx <- NULL
time.trans.idx <- NULL
# '.verbose = TRUE' is used for debugging
# 'when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0' means when central node does not have any neighbour in the mutual info net
local.unrolled.DBN.adj.matrix.list <-
foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
# when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %:%
foreach::foreach(time.trans.idx = 1:num.time.trans, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
when(base::sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %do% # %dopar%
{
centralNodeName <- base::rownames(mi.net.adj.matrix)[centralNodeIdx]
# List names of the central node's neighbours in mi.net.adj.matrix
nbghNames <- base::c()
if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
{
n <- base::nrow(mi.net.adj.matrix)
for (nbrIdx in 1:n)
{
if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
{
nbghNames <- base::rownames(mi.net.adj.matrix)[nbrIdx]
break
}
}
base::rm(n)
}
else if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
{
nbghNames <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
}
local.net.node.names <- base::c(centralNodeName, nbghNames)
# Select subdataset corr. to the nodes in 'local.net.node.names' and time points 'c(time.trans.idx, (time.trans.idx + 1))'
local.input.data.3D <- input.data.discr.3D[base::c(time.trans.idx, (time.trans.idx + 1)), local.net.node.names, ]
#---------------------------------
# Begin: Local Unrolled DBN struct learning
#---------------------------------
## Begin: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# Say, the local nodes are {v1, v2, v3} and the time points are {t1, t2}.
# Then 'local.DBN.input.data.var.names' contains {v1_t1, v2_t1, v3_t1, v1_t2, v2_t2, v3_t2}.
# In that case, 'local.DBN.input.data' will contain columns corr. to elements in 'local.DBN.input.data.var.names' and
# rows corr. to different samples.
## Begin: Generate first row of 'local.DBN.input.data'
sampleIdx <- 1
local.DBN.input.data <- base::matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
ncol = base::dim(local.input.data.3D)[2])
data.to.combine <- base::matrix(local.input.data.3D[2, , sampleIdx], nrow = 1,
ncol = base::dim(local.input.data.3D)[2])
local.DBN.input.data <- base::cbind(local.DBN.input.data, data.to.combine)
base::rm(data.to.combine)
# for (time.pt in 2:dim(local.input.data.3D)[1])
# {
# data.to.combine <- matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# ncol = dim(local.input.data.3D)[2])
# local.DBN.input.data <- cbind(local.DBN.input.data, data.to.combine)
# }
local.DBN.input.data.var.names <- base::c()
for (time.pt in time.trans.idx:(time.trans.idx + 1)) {
for (var.name in base::dimnames(local.input.data.3D)[2]) {
local.DBN.input.data.var.names <- base::c(local.DBN.input.data.var.names,
base::paste(var.name, base::as.character(time.pt), sep = "_t"))
}
}
base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
## End: Generate first row of 'local.DBN.input.data'
## Begin: Generate rest of the rows of 'local.DBN.input.data'
if (base::dim(local.input.data.3D)[3] > 1) # If there are multiple samples per time pt
{
for (sampleIdx in 2:base::dim(local.input.data.3D)[3])
{
sample.to.combine <- base::matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
ncol = base::dim(local.input.data.3D)[2])
data.to.combine <- base::matrix(local.input.data.3D[2, , sampleIdx], nrow = 1,
ncol = base::dim(local.input.data.3D)[2])
sample.to.combine <- base::cbind(sample.to.combine, data.to.combine)
# rm(data.to.combine)
# for (time.pt in 2:dim(local.input.data.3D)[1])
# {
# data.to.combine <- matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# ncol = dim(local.input.data.3D)[2])
# sample.to.combine <- cbind(sample.to.combine, data.to.combine)
# }
# local.DBN.input.data.var.names <- c()
# for (time.pt in 1:dim(local.input.data.3D)[1]) {
# for (var.name in dimnames(local.input.data.3D)[2]) {
# local.DBN.input.data.var.names <- c(local.DBN.input.data.var.names,
# paste(var.name, as.character(time.pt), sep = "_t"))
# }
# }
# colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
local.DBN.input.data <- base::rbind(local.DBN.input.data, sample.to.combine)
}
}
base::rownames(local.DBN.input.data) <- base::as.vector(base::unlist(base::dimnames(local.input.data.3D)[3]))
## End: Generate rest of the rows of 'local.DBN.input.data'
## End: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
# discreteness = rep(TRUE, ncol(local.DBN.input.data)),
# variables = colnames(local.DBN.input.data),
# node.sizes = rep(num.discr.levels, ncol(local.DBN.input.data)),
# starts.from = 0,
# num.time.steps = nrow(local.input.data.3D))
# Begin: Uncomment this section after testing
local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
discreteness = base::rep(TRUE, base::ncol(local.DBN.input.data)),
variables = base::colnames(local.DBN.input.data),
node.sizes = base::rep(num.discr.levels, base::ncol(local.DBN.input.data)),
num.time.steps = base::nrow(local.input.data.3D))
# algo = "mmhc", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# algo = "sm", scoring.func = "BIC", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# scoring.func = 'BIC',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# algo = "mmhc", scoring.func = "BDeu", with layering, initial.network = 'random.chain', seed = 12345
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
layers <- base::c(base::rep(1, base::dim(local.input.data.3D)[2]), base::rep(2, base::dim(local.input.data.3D)[2]))
local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
initial.network = 'random.chain',
seed = 12345,
algo = 'mmhc',
scoring.func = 'BDeu',
num.time.steps =
bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
layering = layers)
# algo = "mmhc", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# # The following four lines must be executed at the same time
# save.plot.to.filename = paste(paste('LocalUnrolledDbn', centralNodeName, sep = '_'), '.jpg', sep = '')
# jpeg(file = paste('LocalUnrolledDbn_', centralNodeName, '.jpg', sep = ''))
# plot(local.unrolled.DBN)
# dev.off()
# Extracting the adjacency matrix of the local DBN
local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
nrow = base::length(local.unrolled.DBN@variables),
ncol = base::length(local.unrolled.DBN@variables),
dimnames = base::c(base::list(local.unrolled.DBN@variables), base::list(local.unrolled.DBN@variables)))
# End: Uncomment this section after testing
# # Begin: This section is for testing
# local.unrolled.DBN.adj.matrix <- matrix(0,
# nrow = ncol(local.DBN.input.data),
# ncol = ncol(local.DBN.input.data),
# dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# # End: This section is for testing
# 'fixed = FALSE' represents that the given pattern is a regular expression.
# local.unrolled.DBN.central.node.indices <- grep(paste('^', centralNodeName, sep = ''),
# colnames(local.unrolled.DBN.adj.matrix),
# fixed = FALSE)
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
local.unrolled.DBN.src.node.names <- local.unrolled.DBN@variables[1:base::dim(local.input.data.3D)[2]]
local.unrolled.DBN.tgt.node.name <- base::paste(centralNodeName, base::as.character((time.trans.idx + 1)), sep = "_t")
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
local.unrolled.DBN.adj.submatrix <- base::matrix(
local.unrolled.DBN.adj.matrix[local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name],
nrow = base::length(local.unrolled.DBN.src.node.names),
ncol = 1,
dimnames = base::c(base::list(local.unrolled.DBN.src.node.names), base::list(local.unrolled.DBN.tgt.node.name)))
# print(local.unrolled.DBN.adj.submatrix)
#---------------------------------
# End: Local Unrolled DBN struct learning
#---------------------------------
# Return value for each 'foreach' iteration
local.unrolled.DBN.adj.submatrix
}
base::rm(time.trans.idx)
base::rm(centralNodeIdx)
# Shut down the cluster
parallel::stopCluster(cl)
# Begin: Unrolled DBN struct learning
# Begin: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
unrolled.DBN.adj.matrix.node.names <- base::c()
for (time.pt in 1:num.timepts) {
for (node.name in base::dimnames(input.data.discr.3D)[2]) {
unrolled.DBN.adj.matrix.node.names <- base::c(unrolled.DBN.adj.matrix.node.names,
base::paste(node.name, base::as.character(time.pt), sep = "_t"))
}
}
unrolled.DBN.adj.matrix <- base::matrix(0, nrow = base::length(unrolled.DBN.adj.matrix.node.names),
ncol = base::length(unrolled.DBN.adj.matrix.node.names),
dimnames = base::c(base::list(unrolled.DBN.adj.matrix.node.names),
base::list(unrolled.DBN.adj.matrix.node.names)))
# End: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
# 'local.unrolled.DBN.adj.matrix.list' is a list of lists of matrices.
# The outer list contains 'local.unrolled.DBN.adj.matrix.list' inner lists.
# Each inner list contains 'num.time.trans' matrices.
# length(local.unrolled.DBN.adj.matrix.list) = num.nodes
for (outer.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list))
{
# length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = num.time.trans
for (inner.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]))
{
# 'submatrix.to.combine'
submatrix.to.combine <- local.unrolled.DBN.adj.matrix.list[[outer.list.idx]][[inner.list.idx]]
unrolled.DBN.adj.matrix[base::rownames(submatrix.to.combine), base::colnames(submatrix.to.combine)] <- submatrix.to.combine
}
}
# End: Unrolled DBN struct learning
# Return the final directed net of the current iteration
return (unrolled.DBN.adj.matrix)
}
#' Unrolled DBN structure learning with Markov Order 1
#'
#' Candidate parents: The target node itself and its CLR net neighbours at immediately previous time pt.
#'
#' @param input.data.discr.3D Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
#' @param mi.net.adj.matrix Adjacency matrix of the mutual information network. Rownames and colnames should be node names.
#' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized for performing the DBN structure learning.
#' @param num.nodes number of nodes
#' @param num.timepts number of timepoints
#' @param max.fanin maximum incoming edges in the graph
#' @param output.dirname output directory to store files
#'
#' @return Unrolled DBN adjacency matrix
#'
#' @keywords internal
#' @noRd
learnDbnStructMo1Layer3dParDeg1 <- function(input.data.discr.3D, mi.net.adj.matrix, num.discr.levels, num.nodes, num.timepts, max.fanin, output.dirname) {
if(!base::is.array(input.data.discr.3D))
{
base::stop("Error in learnDbnStructMo1Layer3dParDeg1 input.data.discr.3D is not an array")
}
if(!base::is.matrix(mi.net.adj.matrix))
{
base::stop("Error in learnDbnStructMo1Layer3dParDeg1 mi.net.adj.matrix is not a matrix")
}
num.time.trans <- (num.timepts - 1)
## Start and register a parallel backend for parallel computing
## 10 cores to be used for grni server
# no_cores <- min(10, num.nodes, (parallel::detectCores() - 1))
# cl <- parallel::makeCluster(no_cores, outfile = paste(getwd(), 'asset/outfile.txt', sep = '/' ))
# doParallel::registerDoParallel(cl)
# '.verbose = TRUE' is used for debugging
# 'when(sum(mi.net.adj.matrix[, centralNodeIdx]) > 0' means when central node has at least one neighbour in the mutual info net
# Use %do% amd %dopar% for serial and parallel computing, resp.
centralNodeIdx <- NULL
time.trans.idx <- NULL
local.unrolled.DBN.adj.matrix.list <-
foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
# when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %:%
foreach::foreach(time.trans.idx = 1:num.time.trans, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
when(base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 0) %do%
{
centralNodeName <- base::rownames(mi.net.adj.matrix)[centralNodeIdx]
# List names of the central node's neighbours in mi.net.adj.matrix
nbghNames <- base::c()
if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
{
for (nbrIdx in 1:base::nrow(mi.net.adj.matrix))
{
if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
{
nbghNames <- base::rownames(mi.net.adj.matrix)[nbrIdx]
break
}
}
}
else if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
{
nbghNames <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
}
local.net.node.names <- base::c(centralNodeName, nbghNames)
#---------------------------------
# Begin: Local Unrolled DBN struct learning
#---------------------------------
## Begin: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# Say, the central node is v1 and the local nodes are {v1, v2, v3}, and the time points are {t1, t2}.
# Then 'local.DBN.input.data.var.names' contains {v1_t1, v2_t1, v3_t1, v1_t2}.
# In that case, 'local.DBN.input.data' will contain columns corr. to elements in 'local.DBN.input.data.var.names' and
# rows corr. to different samples.
num.samples <- base::dim(input.data.discr.3D)[3]
local.DBN.input.data <- base::matrix(0, nrow = num.samples, ncol = (base::length(local.net.node.names) + 1))
local.unrolled.DBN.src.node.names <- base::c()
for (var.name in local.net.node.names) {
local.unrolled.DBN.src.node.names <- base::c(local.unrolled.DBN.src.node.names,
base::paste(var.name, base::as.character(time.trans.idx), sep = "_t"))
}
local.unrolled.DBN.tgt.node.name <- base::paste(centralNodeName, base::as.character(time.trans.idx + 1), sep = "_t")
local.DBN.input.data.var.names <- base::c(local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name)
base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
for (node.name in local.net.node.names)
{
col.name <- base::paste(node.name, base::as.character(time.trans.idx), sep = "_t")
local.DBN.input.data[, col.name] <- input.data.discr.3D[time.trans.idx, node.name, ]
}
col.name <- base::paste(centralNodeName, base::as.character(time.trans.idx + 1), sep = "_t")
local.DBN.input.data[, col.name] <- input.data.discr.3D[(time.trans.idx + 1), centralNodeName, ]
base::rm(col.name)
## End: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
discreteness = base::rep(TRUE, base::ncol(local.DBN.input.data)),
variables = base::colnames(local.DBN.input.data),
node.sizes = base::rep(num.discr.levels, base::ncol(local.DBN.input.data)))
# algo = "mmhc", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, (ncol(local.DBN.input.data) - 1)), 2)
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# rm(layers)
## algo = "sm", scoring.func = "BIC", with layering
## There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
## The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
## A node with layer idx j can have parents from layer idx i such that i =< j.
layers <- base::c(base::rep(1, (base::ncol(local.DBN.input.data) - 1)), 2)
local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
algo = 'sm',
scoring.func = 'BIC',
layering = layers)
base::rm(layers)
## algo = "sm", scoring.func = "BIC", with layering, with max.fanin
## There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
## The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
## A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, (ncol(local.DBN.input.data) - 1)), 2)
# local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# scoring.func = 'BIC',
# max.fanin = max.fanin,
# layering = layers)
# rm(layers)
# algo = "mmhc", scoring.func = "BDeu", with layering, initial.network = 'random.chain', seed = 12345
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# initial.network = 'random.chain',
# seed = 12345,
# algo = 'mmhc',
# scoring.func = 'BDeu',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# algo = "mmhc", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# # The following four lines must be executed at the same time
# save.plot.to.filename = paste(paste('LocalUnrolledDbn', centralNodeName, sep = '_'), '.jpg', sep = '')
# jpeg(file = paste('LocalUnrolledDbn_', centralNodeName, '.jpg', sep = ''))
# plot(local.unrolled.DBN)
# dev.off()
# Extracting the adjacency matrix of the local DBN
local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
nrow = base::length(local.unrolled.DBN@variables),
ncol = base::length(local.unrolled.DBN@variables),
dimnames = base::c(base::list(local.unrolled.DBN@variables), base::list(local.unrolled.DBN@variables)))
# This for loop checks whether parents are learnt only for 'local.unrolled.DBN.tgt.node.name'. If
# so, then nothing is printed. Otherwise, prints the column(s) corr. to the undesired tgt node(s).
for (col.idx in 1:(base::ncol(local.unrolled.DBN.adj.matrix) - 1))
{
if (base::sum(local.unrolled.DBN.adj.matrix[, col.idx]) > 0)
{
base::print('Erroneous column')
base::print(local.unrolled.DBN.adj.matrix[, col.idx])
}
}
# End: Uncomment this section after testing
# # Begin: This section is for testing
# local.unrolled.DBN.adj.matrix <- matrix(0,
# nrow = ncol(local.DBN.input.data),
# ncol = ncol(local.DBN.input.data),
# dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# # End: This section is for testing
# 'fixed = FALSE' represents that the given pattern is a regular expression.
# local.unrolled.DBN.central.node.indices <- grep(paste('^', centralNodeName, sep = ''),
# colnames(local.unrolled.DBN.adj.matrix),
# fixed = FALSE)
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
local.unrolled.DBN.adj.submatrix <- base::matrix(
local.unrolled.DBN.adj.matrix[local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name],
nrow = base::length(local.unrolled.DBN.src.node.names),
ncol = 1,
dimnames = base::c(base::list(local.unrolled.DBN.src.node.names), base::list(local.unrolled.DBN.tgt.node.name)))
# print(local.unrolled.DBN.adj.submatrix)
#---------------------------------
# End: Local Unrolled DBN struct learning
#---------------------------------
# Return value for each 'foreach' iteration
local.unrolled.DBN.adj.submatrix
}
base::rm(centralNodeIdx)
base::rm(time.trans.idx)
##print('End of foreach loops')
## Shut down the cluster
# parallel::stopCluster(cl)
# Begin: Unrolled DBN struct learning
# Begin: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
unrolled.DBN.adj.matrix.node.names <- base::c()
for (time.pt in 1:num.timepts) {
for (node.name in base::dimnames(input.data.discr.3D)[2]) {
unrolled.DBN.adj.matrix.node.names <- base::c(unrolled.DBN.adj.matrix.node.names,
base::paste(node.name, base::as.character(time.pt), sep = "_t"))
}
}
unrolled.DBN.adj.matrix <- base::matrix(0, nrow = base::length(unrolled.DBN.adj.matrix.node.names),
ncol = base::length(unrolled.DBN.adj.matrix.node.names),
dimnames = base::c(base::list(unrolled.DBN.adj.matrix.node.names),
base::list(unrolled.DBN.adj.matrix.node.names)))
# End: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
## Debugging
# writeLines('length(local.unrolled.DBN.adj.matrix.list) = ', length(local.unrolled.DBN.adj.matrix.list), '\n')
# print(length(local.unrolled.DBN.adj.matrix.list))
# 'local.unrolled.DBN.adj.matrix.list' is a list of lists of matrices.
# The outer list contains 'local.unrolled.DBN.adj.matrix.list' inner lists.
# Each inner list contains 'num.time.trans' matrices.
# length(local.unrolled.DBN.adj.matrix.list) = num.nodes
for (outer.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list))
{
## Debugging
# writeLines('outer.list.idx = ', outer.list.idx, '\n')
# print(outer.list.idx)
# writeLines('length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = ',
# length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]), '\n')
# print(length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]))
# print(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]])
# if 'local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]' is not an empty list
if (base::length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) > 0)
{
# length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = num.time.trans
for (inner.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]))
{
## Debugging
# writeLines('inner.list.idx = ', inner.list.idx, '\n')
# print(inner.list.idx)
# 'submatrix.to.combine'
submatrix.to.combine <- local.unrolled.DBN.adj.matrix.list[[outer.list.idx]][[inner.list.idx]]
unrolled.DBN.adj.matrix[base::rownames(submatrix.to.combine), base::colnames(submatrix.to.combine)] <- submatrix.to.combine
}
}
}
# End: Unrolled DBN struct learning
# Return the final directed net of the current iteration
return (unrolled.DBN.adj.matrix)
}
#' Unrolled DBN structure learning with Markov Order 1
#'
#' This function is the newer version of learnDbnStructMo1Layer3dParDeg1().
#' The only difference is in the size of unrolled DBN adjacency matrix. In earlier version,
#' the size is ((V \ times T) \ times (v \ times T)) where V = number of nodes and
#' T = number of time points. But the size is too large when V is very large. For
#' e.g., when V = 4028, the function can not execute successfully even with 32 GB main
#' memory in grni server. In this version, the size is reduced by storing the unrolled DBN
#' in an adjacency list of length (T - 1). The t^{th} element in the list is
#' the predicted network adjacency matrix at the t^{th} time interval. Therefore, each list element
#' is a binary matrix of dimension (V \ times V). Hence, the total size of the unrolled DBN
#' adjacency list is ((T - 1) \ times (V \ times V)).
#' Candidate parents: The target node itself and its CLR net neighbours at immediately previous time pt.
#'
#' @param input.data.discr.3D Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
#' @param mi.net.adj.matrix Adjacency matrix of the mutual information network. Rownames and colnames should be node names.
#' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized for performing the DBN structure learning.
#' @param num.nodes number of nodes
#' @param num.timepts number of timepoints
#' @param max.fanin maximum incoming edges in the graph
#' @param node.names name of the nodes
#' @param clr.algo clr algorithm to use
#'
#' @return Unrolled DBN adjacency matrix
#'
#' @keywords internal
#' @noRd
learnDbnStructMo1Layer3dParDeg1_v2 <- function(input.data.discr.3D, mi.net.adj.matrix,
num.discr.levels, num.nodes, num.timepts, max.fanin,
node.names, clr.algo) {
if(!base::is.array(input.data.discr.3D))
{
base::stop("Error in learnDbnStructMo1Layer3dParDeg1_v2 input.data.discr.3D is not an array")
}
if(!base::is.matrix(mi.net.adj.matrix))
{
base::stop("Error in learnDbnStructMo1Layer3dParDeg1_v2 mi.net.adj.matrix is not a matrix")
}
num.time.trans <- (num.timepts - 1)
## Start and register a parallel backend for parallel computing
## 10 cores to be used for grni server
# no_cores <- min(10, num.nodes, (parallel::detectCores() - 1))
# cl <- parallel::makeCluster(no_cores, outfile = paste(getwd(), 'asset/outfile.txt', sep = '/' ))
# doParallel::registerDoParallel(cl)
# '.verbose = TRUE' is used for debugging
# 'when(sum(mi.net.adj.matrix[, centralNodeIdx]) > 0' means when central node has at least one neighbour in the mutual info net
# Use %do% amd %dopar% for serial and parallel computing, resp.
centralNodeIdx <-NULL
time.trans.idx <-NULL
local.unrolled.DBN.adj.matrix.list <-
foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
# when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %:%
foreach::foreach(time.trans.idx = 1:num.time.trans, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
when(base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 0) %do%
{
centralNodeName <- base::rownames(mi.net.adj.matrix)[centralNodeIdx]
# List names of the central node's neighbours in mi.net.adj.matrix
nbghNames <-base::c()
if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
{
for (nbrIdx in 1:base::nrow(mi.net.adj.matrix))
{
if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
{
nbghNames <- base::rownames(mi.net.adj.matrix)[nbrIdx]
break
}
}
}
else if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
{
nbghNames <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
}
local.net.node.names <- base::c()
if (clr.algo == 'CLR2.1') {
local.net.node.names <- nbghNames
} else if ((clr.algo == 'CLR') | (clr.algo == 'CLR2')) {
local.net.node.names <- base::c(centralNodeName, nbghNames)
}
#---------------------------------
# Begin: Local Unrolled DBN struct learning
#---------------------------------
## Begin: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# Say, the central node is v1 and the local nodes are {v1, v2, v3}, and the time points are {t1, t2}.
# Then 'local.DBN.input.data.var.names' contains {v1_t1, v2_t1, v3_t1, v1_t2}.
# In that case, 'local.DBN.input.data' will contain columns corr. to elements in 'local.DBN.input.data.var.names' and
# rows corr. to different samples.
num.samples <- base::dim(input.data.discr.3D)[3]
local.DBN.input.data <- base::matrix(0, nrow = num.samples, ncol = (base::length(local.net.node.names) + 1))
local.unrolled.DBN.src.node.names <- base::c()
for (var.name in local.net.node.names) {
local.unrolled.DBN.src.node.names <- base::c(local.unrolled.DBN.src.node.names,
base::paste(var.name, base::as.character(time.trans.idx), sep = "_t"))
}
base::rm(var.name)
local.unrolled.DBN.tgt.node.name <- base::paste(centralNodeName, base::as.character(time.trans.idx + 1), sep = "_t")
local.DBN.input.data.var.names <- base::c(local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name)
base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
for (node.name in local.net.node.names) {
col.name <- base::paste(node.name, base::as.character(time.trans.idx), sep = "_t")
local.DBN.input.data[, col.name] <- input.data.discr.3D[time.trans.idx, node.name, ]
}
base::rm(node.name)
col.name <- base::paste(centralNodeName, base::as.character(time.trans.idx + 1), sep = "_t")
local.DBN.input.data[, col.name] <- input.data.discr.3D[(time.trans.idx + 1), centralNodeName, ]
base::rm(col.name)
## End: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# save(local.DBN.input.data, file = 'local.DBN.input.data.RData')
local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
discreteness = base::rep(TRUE, ncol(local.DBN.input.data)),
variables = base::colnames(local.DBN.input.data),
node.sizes = base::rep(num.discr.levels,
base::ncol(local.DBN.input.data)))
# algo = "mmhc", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, (ncol(local.DBN.input.data) - 1)), 2)
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# rm(layers)
## algo = "sm", scoring.func = "BIC", with layering
## There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
## The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
## A node with layer idx j can have parents from layer idx i such that i =< j.
layers <- base::c(base::rep(1, (base::ncol(local.DBN.input.data) - 1)), 2)
local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
algo = 'sm',
scoring.func = 'BIC',
layering = layers)
base::rm(layers)
## algo = "sm", scoring.func = "BIC", with layering, with max.fanin
## There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
## The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
## A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, (ncol(local.DBN.input.data) - 1)), 2)
# local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# scoring.func = 'BIC',
# max.fanin = max.fanin,
# layering = layers)
# rm(layers)
# algo = "mmhc", scoring.func = "BDeu", with layering, initial.network = 'random.chain', seed = 12345
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# initial.network = 'random.chain',
# seed = 12345,
# algo = 'mmhc',
# scoring.func = 'BDeu',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# algo = "mmhc", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# # The following four lines must be executed at the same time
# save.plot.to.filename = paste(paste('LocalUnrolledDbn', centralNodeName, sep = '_'), '.jpg', sep = '')
# jpeg(file = paste('LocalUnrolledDbn_', centralNodeName, '.jpg', sep = ''))
# plot(local.unrolled.DBN)
# dev.off()
# Extracting the adjacency matrix of the local DBN
local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
nrow = base::length(local.unrolled.DBN@variables),
ncol = base::length(local.unrolled.DBN@variables),
dimnames = base::c(base::list(local.unrolled.DBN@variables), base::list(local.unrolled.DBN@variables)))
# This for loop checks whether parents are learnt only for 'local.unrolled.DBN.tgt.node.name'. If
# so, then nothing is printed. Otherwise, prints the column(s) corr. to the undesired tgt node(s).
for (col.idx in 1:(base::ncol(local.unrolled.DBN.adj.matrix) - 1))
{
if (base::sum(local.unrolled.DBN.adj.matrix[, col.idx]) > 0)
{
base::print('Erroneous column')
base::print(local.unrolled.DBN.adj.matrix[, col.idx])
}
}
# End: Uncomment this section after testing
# # Begin: This section is for testing
# local.unrolled.DBN.adj.matrix <- matrix(0,
# nrow = ncol(local.DBN.input.data),
# ncol = ncol(local.DBN.input.data),
# dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# # End: This section is for testing
# 'fixed = FALSE' represents that the given pattern is a regular expression.
# local.unrolled.DBN.central.node.indices <- grep(paste('^', centralNodeName, sep = ''),
# colnames(local.unrolled.DBN.adj.matrix),
# fixed = FALSE)
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
## Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
local.unrolled.DBN.adj.submatrix <- base::matrix(
local.unrolled.DBN.adj.matrix[local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name],
nrow = base::length(local.unrolled.DBN.src.node.names),
ncol = 1,
dimnames = base::c(base::list(local.unrolled.DBN.src.node.names), base::list(local.unrolled.DBN.tgt.node.name)))
# print(local.unrolled.DBN.adj.submatrix)
#---------------------------------
# End: Local Unrolled DBN struct learning
#---------------------------------
# Return value for each 'foreach' iteration
local.unrolled.DBN.adj.submatrix
}
base::rm(centralNodeIdx)
base::rm(time.trans.idx)
##print('End of foreach loops')
# save(local.unrolled.DBN.adj.matrix.list, file = paste(getwd(), 'asset/local.unrolled.DBN.adj.matrix.list.RData', sep = '/'))
## Shut down the cluster
# parallel::stopCluster(cl)
# Begin: Unrolled DBN struct learning
# Initialize the unrolled DBN adjacency list.
## It is a list of length (T - 1) where T = number of time points. The t^{th} element in the list is
## the predicted network adjacency matrix at the t^{th} time interval. Therefore, each list element
## is a binary matrix of dimension (V \times V) where V = number of nodes. Hence, the total size
## of the unrolled DBN adjacency list is ((T - 1) \times (V \times V)).
## Each adjacency matrix is initialized with a zero matrix of dimension (V \times V). Its rows
## corr. to soruce nodes and the columns corr. to target nodes.
unrolled.DBN.adj.matrix.list <- base::list()
for (list.idx in 1:num.time.trans)
{
unrolled.DBN.adj.matrix.list[[list.idx]] <- base::matrix(0, nrow = num.nodes,
ncol = num.nodes,
dimnames =
base::c(base::list(node.names),
base::list(node.names)))
}
base::rm(list.idx)
## 'local.unrolled.DBN.adj.matrix.list' is a list of lists of matrices.
## The outer list contains length(local.unrolled.DBN.adj.matrix.list) number of inner lists.
## The length is \ge 1 and \le num.nodes. It is < num.nodes when there exists some node
## without any neighbour in the mutual info net.
##
## Each inner list contains 'num.time.trans' matrices.
for (outer.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list))
{
## Debugging
# writeLines('outer.list.idx = ', outer.list.idx, '\n')
# print(outer.list.idx)
# writeLines('length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = ',
# length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]), '\n')
# print(length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]))
# print(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]])
## if 'local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]' is not an empty list
if (base::length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) > 0)
{
## length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = num.time.trans
## for any value of outer.list.idx
for (inner.list.idx in 1:num.time.trans)
{
## Debugging
# writeLines('inner.list.idx = ', inner.list.idx, '\n')
# print(inner.list.idx)
## 'submatrix.to.combine'
submatrix.to.combine <- local.unrolled.DBN.adj.matrix.list[[outer.list.idx]][[inner.list.idx]]
## Begin: remove the timestamps from the row and col names of submatrix.to.combine.
## E.g., a row or col name 'G2_t34' will be converted to just 'G2'.
submatrix.to.combine.rownames <- base::c()
for (row.idx in 1:base::nrow(submatrix.to.combine))
{
old.rowname <- base::rownames(submatrix.to.combine)[row.idx]
## Substitute '_t[0-9]+' pattern with empty string in old.rowname
new.rowname <- base::sub('_t[0-9]+', '', old.rowname)
submatrix.to.combine.rownames <- base::c(submatrix.to.combine.rownames, new.rowname)
}
base::rm(row.idx)
base::rownames(submatrix.to.combine) <- submatrix.to.combine.rownames
base::rm(submatrix.to.combine.rownames)
submatrix.to.combine.colnames <- base::c()
for (col.idx in 1:base::ncol(submatrix.to.combine))
{
old.colname <- base::colnames(submatrix.to.combine)[col.idx]
## Replace '_t[0-9]+' pattern with empty string in old.rowname
new.colname <- base::sub('_t[0-9]+', '', old.colname)
submatrix.to.combine.colnames <-base::c(submatrix.to.combine.colnames, new.colname)
}
base::rm(col.idx)
base::colnames(submatrix.to.combine) <- submatrix.to.combine.colnames
base::rm(submatrix.to.combine.colnames)
## End: remove the timestamps from the row and col names of submatrix.to.combine.
# print('submatrix.to.combine')
# print(submatrix.to.combine)
#
# print('inner.list.idx')
# print(inner.list.idx)
# print('unrolled.DBN.adj.matrix.list[[inner.list.idx]]')
# print(unrolled.DBN.adj.matrix.list[[inner.list.idx]])
unrolled.DBN.adj.matrix.list[[inner.list.idx]][base::rownames(submatrix.to.combine),
base::colnames(submatrix.to.combine)] <- submatrix.to.combine
}
}
}
base::rm(outer.list.idx)
# End: Unrolled DBN struct learning
return (unrolled.DBN.adj.matrix.list)
}
#' Learns DBN structure of Markov order 1 where candidate parents are selected using the CLR3 algo
#'
#' It is a serial algorithmic implementation.
#'
#' @param input.data.discr.3D Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
#' @param mi.net.adj.matrix.list.filename List of filenames (relative/absolute path) of adjacency matrices of the mutual information network.
#' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized for performing the DBN structure learning.
#' @param num.nodes number of nodes
#' @param num.timepts number of timepoints
#' @param max.fanin maximum incoming edges in the graph
#' @param node.names name of the nodes
#' @param unrolled.DBN.adj.matrix.list List of unrolled DBN adjacency matrices
#'
#' @return Unrolled DBN adjacency matrix
#'
#' @keywords internal
#' @noRd
LearnDbnStructMo1Clr3Ser <- function(input.data.discr.3D, mi.net.adj.matrix.list.filename,
num.discr.levels, num.nodes, num.timepts, max.fanin,
node.names, unrolled.DBN.adj.matrix.list) {
if(!base::is.array(input.data.discr.3D))
{
base::stop("Error in LearnDbnStructMo1Clr3Ser input.data.discr.3D is not an array")
}
if(!base::is.list(mi.net.adj.matrix.list.filename))
{
base::stop("Error in LearnDbnStructMo1Clr3Ser mi.net.adj.matrix.list.filename is not a list")
}
if(!base::is.list(unrolled.DBN.adj.matrix.list))
{
base::stop("Error in LearnDbnStructMo1Clr3Ser unrolled.DBN.adj.matrix.list is not a list")
}
## Here, each 'time.pt.idx' represents time interval
## ('time.pt.idx', ('time.pt.idx' + 1))
for (time.pt.idx in 1:(num.timepts - 1)) {
## Load the adjacency matrices of time-interval-specific
## CLR nets from the given file
mi.net.adj.matrix.list <- NULL
base::load(mi.net.adj.matrix.list.filename)
## Extract the CLR net specific to the current time interval
mi.net.adj.matrix <- mi.net.adj.matrix.list[[time.pt.idx]]
base::rm(mi.net.adj.matrix.list.filename)
for (tgt.node.idx in 1:num.nodes) {
## If the current target node has no candidate parent,
## then skip to the next target node
if (base::sum(mi.net.adj.matrix[, tgt.node.idx]) == 0) {
next
}
## E.g., 'v1_t2'
## Do not use rownames
tgt.node.name <- base::colnames(mi.net.adj.matrix)[tgt.node.idx]
## E.g., 'v1'
tgt.node.base.name <- node.names[tgt.node.idx]
tgt.node.name.as.src <- base::paste(tgt.node.base.name, base::as.character(time.pt.idx), sep = "_t")
## List names of the target node's neighbours in 'mi.net.adj.matrix'
ngbr.names <- base::c()
if (base::sum(mi.net.adj.matrix[, tgt.node.idx]) == 1) {
## Just one neighbour
for (ngbr.idx in 1:base::nrow(mi.net.adj.matrix)) {
if (mi.net.adj.matrix[ngbr.idx, tgt.node.idx] == 1) {
## Do not use colnames
ngbr.names <- base::rownames(mi.net.adj.matrix)[ngbr.idx]
break
}
}
base::rm(ngbr.idx)
} else if (base::sum(mi.net.adj.matrix[, tgt.node.idx]) > 1) {
## Multiple neighbours
## Do not use colnames
ngbr.names <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, tgt.node.idx] == 1),])
}
local.net.node.names <- base::c(tgt.node.name.as.src, ngbr.names)
##------------------------------------------------------------
## Begin: Local Unrolled DBN struct learning
##------------------------------------------------------------
num.samples.per.timept <- base::dim(input.data.discr.3D)[3]
local.DBN.input.data <- base::matrix(0, nrow = num.samples.per.timept,
ncol = (base::length(local.net.node.names) + 1))
local.unrolled.DBN.src.node.names <- base::c()
local.unrolled.DBN.src.node.names <- base::c(ngbr.names, tgt.node.name.as.src)
local.unrolled.DBN.src.node.names <- base::sort(local.unrolled.DBN.src.node.names)
base::rm(tgt.node.name.as.src)
local.unrolled.DBN.tgt.node.name <- tgt.node.name
local.DBN.input.data.var.names <- base::c(local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name)
base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
base::rm(local.DBN.input.data.var.names)
for (node.name in local.net.node.names) {
## Substitute '_t[0-9]+' pattern with empty string in 'node.name'.
## If 'node.name' = 'v1_t1', then
## 'node.base.name' = 'v1'.
node.base.name <- base::sub('_t[0-9]+', '', node.name)
local.DBN.input.data[, node.name] <- input.data.discr.3D[time.pt.idx, node.base.name, ]
}
base::rm(node.name)
local.DBN.input.data[, tgt.node.name] <- input.data.discr.3D[(time.pt.idx + 1), tgt.node.base.name, ]
local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
discreteness = base::rep(TRUE, base::ncol(local.DBN.input.data)),
variables = base::colnames(local.DBN.input.data),
node.sizes = base::rep(num.discr.levels,
base::ncol(local.DBN.input.data)))
layers <- base::c(base::rep(1, (base::ncol(local.DBN.input.data) - 1)), 2)
local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
algo = 'sm',
scoring.func = 'BIC',
layering = layers)
base::rm(layers)
## Extracting the adjacency matrix of the local DBN
local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
nrow = base::length(local.unrolled.DBN@variables),
ncol = base::length(local.unrolled.DBN@variables),
dimnames = base::c(base::list(local.unrolled.DBN@variables),
base::list(local.unrolled.DBN@variables)))
# This for loop checks whether parents are learnt only for 'local.unrolled.DBN.tgt.node.name'. If
# so, then nothing is printed. Otherwise, prints the column(s) corr. to the undesired tgt node(s).
for (col.idx in 1:(base::ncol(local.unrolled.DBN.adj.matrix) - 1)) {
if (base::sum(local.unrolled.DBN.adj.matrix[, col.idx]) > 0) {
base::print('Erroneous column')
base::print(local.unrolled.DBN.adj.matrix[, col.idx])
}
}
base::rm(col.idx)
## Assuming there are > 1 time points.
## Otherwise, 'submatrix.to.combine' would become a vector.
submatrix.to.combine <- base::matrix(
local.unrolled.DBN.adj.matrix[local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name],
nrow = base::length(local.unrolled.DBN.src.node.names),
ncol = 1,
dimnames = base::c(base::list(local.unrolled.DBN.src.node.names), base::list(local.unrolled.DBN.tgt.node.name)))
##------------------------------------------------------------
## End: Local Unrolled DBN struct learning
##------------------------------------------------------------
## Begin: remove the timestamps from the row names of
## 'submatrix.to.combine'.
## E.g., a row name 'G2_t34' will be converted to just 'G2'.
submatrix.to.combine.rownames <- base::c()
for (row.idx in 1:base::nrow(submatrix.to.combine)) {
old.rowname <- base::rownames(submatrix.to.combine)[row.idx]
## Substitute '_t[0-9]+' pattern with empty string in old.rowname
new.rowname <- base::sub('_t[0-9]+', '', old.rowname)
submatrix.to.combine.rownames <- base::c(submatrix.to.combine.rownames, new.rowname)
}
base::rm(row.idx)
base::rownames(submatrix.to.combine) <- submatrix.to.combine.rownames
base::rm(submatrix.to.combine.rownames)
## End: remove the timestamps from the row names of submatrix.to.combine.
unrolled.DBN.adj.matrix.list[[time.pt.idx]][rownames(submatrix.to.combine),
tgt.node.base.name] <- submatrix.to.combine
}
base::rm(tgt.node.idx)
}
base::rm(time.pt.idx)
return(unrolled.DBN.adj.matrix.list)
}
Try the TGS package in your browser
Any scripts or data that you put into this service are public.
TGS documentation built on July 1, 2020, 10:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions LearnDbnStructMo1Clr3Ser learnDbnStructMo1Layer3dParDeg1_v2 learnDbnStructMo1Layer3dParDeg1 learnDbnStructLayer3dParDeg1
# #Goal: Infer Dynamic Bayesian Network (DBN)
#
# #---------------------------------
# # Begin: Step 2: Decomposing the network
# #---------------------------------
#
#
# #' Unrolled DBN structure learning with all possible Markov Orders.
# #'
# #' Unrolled DBN structure learning with all possible Markov Orders.
# #' Candidate parents: The target node itself and its CLR net neighbours at any previous and current time pt.
# #'
# #' @import bnstruct
# #' @import ggm
# #' @import foreach
# #' @import doParallel
# #'
# #' @param input.data.discr.3D 3D input data Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
# #' @param mi.net.adj.matrix Adjacency matrix of the mutual information network. Rownames and colnames should be node names.
# #' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized
# #' @param num.nodes number of nodes
# #' @param num.timepts number of timepoints
# #'
# #' @return Unrolled DBN adjacency matrix
# #'
# #' @export
# learnDbnStruct3dParDeg1 <- function(input.data.discr.3D, mi.net.adj.matrix, num.discr.levels, num.nodes, num.timepts)
# {
# # load('dream3.yeast1.size10.trajectory.3D.Rdata')
# # input.data.3D <- dream3.yeast1.size10.trajectory.3D.data
# # load('mi.net.adj.matrix.Rdata')
#
# # In 'di.net.adj.matrix', rows are src nodes and cols are tgt nodes
# di.net.adj.matrix <- mi.net.adj.matrix
# di.net.adj.matrix[1:base::nrow(di.net.adj.matrix), 1:base::ncol(di.net.adj.matrix)] <- 0
#
# no_cores <- num.nodes
# cl <- parallel::makeCluster(no_cores)
# doParallel::registerDoParallel(cl)
#
# #base::print('going inside foreach')
#
# # For each central node, learn local DBN struct in parallel.
# # Package names, specified by '.packages' must be copied to each worker core.
# # Current environmental variables that are referenced inside the foreach loop are copied
# # to each worker automatically.
# # local.unrolled.DBN.adj.matrix.list <- foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = 'bnstruct') %dopar%
# # local.unrolled.DBN.adj.matrix.list <- as.list(1:10)
#
# # local.unrolled.DBN.adj.matrix.list <- foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = 'bnstruct') %:% when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %dopar%
# # {
# # print('now inside foreach 1')
# #
# # centralNodeName <- rownames(mi.net.adj.matrix)[centralNodeIdx]
# #
# # # List names of the central node's neighbours in mi.net.adj.matrix
# # nbghNames <- c()
# # if (sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
# # {
# # for (nbrIdx in 1:n)
# # {
# # if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
# # {
# # nbghNames <- rownames(mi.net.adj.matrix)[nbrIdx]
# # break
# # }
# # }
# # }
# # else if (sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
# # {
# # nbghNames <- rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
# # }
# #
# # local.net.node.names <- c(centralNodeName, nbghNames)
# #
# # local.input.data.3D <- input.data.discr.3D[, local.net.node.names, ]
# #
# # #---------------------------------
# # # Begin: Local Unrolled DBN struct learning
# # #---------------------------------
# #
# # sampleIdx <- 1
# #
# # local.DBN.input.data <- matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
# # ncol = dim(local.input.data.3D)[2])
# #
# # for (time.pt in 2:dim(local.input.data.3D)[1])
# # {
# # data.to.combine <- matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# # ncol = dim(local.input.data.3D)[2])
# # local.DBN.input.data <- cbind(local.DBN.input.data, data.to.combine)
# # }
# #
# # local.DBN.input.data.var.names <- c()
# # for (time.pt in 1:dim(local.input.data.3D)[1]) {
# # for (var.name in dimnames(local.input.data.3D)[2]) {
# # local.DBN.input.data.var.names <- c(local.DBN.input.data.var.names,
# # paste(var.name, as.character(time.pt), sep = "_t"))
# # }
# # }
# # colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
# #
# # if (dim(local.input.data.3D)[3] > 1) # If there are multiple samples per time pt
# # {
# # for (sampleIdx in 2:dim(local.input.data.3D)[3])
# # {
# # sample.to.combine <- matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
# # ncol = dim(local.input.data.3D)[2])
# #
# # for (time.pt in 2:dim(local.input.data.3D)[1])
# # {
# # data.to.combine <- matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# # ncol = dim(local.input.data.3D)[2])
# # sample.to.combine <- cbind(sample.to.combine, data.to.combine)
# # }
# #
# # local.DBN.input.data <- rbind(local.DBN.input.data, sample.to.combine)
# # }
# # }
# #
# # rownames(local.DBN.input.data) <- as.vector(unlist(dimnames(local.input.data.3D)[3]))
# #
# # local.unrolled.DBN.adj.matrix <- matrix(0, nrow = ncol(local.DBN.input.data), ncol = ncol(local.DBN.input.data),
# # dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# #
# # # 'fixed = FALSE' represents that the given pattern is a regular expression.
# # local.unrolled.DBN.central.node.indices <- grep(paste('^', centralNodeName, sep = ''),
# # colnames(local.unrolled.DBN.adj.matrix),
# # fixed = FALSE)
# #
# # # Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# # local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
# #
# # #---------------------------------
# # # End: Local Unrolled DBN struct learning
# # #---------------------------------
# #
# # local.unrolled.DBN.adj.submatrix
# # }
# #
# # print(local.unrolled.DBN.adj.matrix.list)
#
# local.unrolled.DBN.adj.matrix.list <- foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = 'bnstruct') %:% when(base::sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %dopar%
# {
#
# # print('now inside foreach')
#
# # centralNodeIdx <- 1
#
# centralNodeName <- base::rownames(mi.net.adj.matrix)[centralNodeIdx]
#
# # if central node does not have any neighbour in mutual information net
# if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 0)
# {
# next
# }
#
# # List names of the central node's neighbours in mi.net.adj.matrix
# nbghNames <- base::c()
# if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
# {
# for (nbrIdx in 1:n)
# {
# if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
# {
# nbghNames <- base::rownames(mi.net.adj.matrix)[nbrIdx]
# break
# }
# }
# }
# else if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
# {
# nbghNames <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
# }
#
# local.net.node.names <- base::c(centralNodeName, nbghNames)
#
# local.input.data.3D <- input.data.discr.3D[, local.net.node.names, ]
#
#
# # #---------------------------------
# # # Begin: Local BN struct learning
# # #---------------------------------
# # local.BN.input.data <- bnstruct::BNDataset(local.discr.data,
# # variables = colnames(local.discr.data),
# # discreteness = rep(TRUE, ncol(local.discr.data)),
# # node.sizes = rep(2, ncol(local.discr.data)),
# # starts.from = 0)
# #
# # # Default params: scoring.func = "BDeu"
# # local.BN <- bnstruct::learn.network(local.discr.data, algo = 'sm')
# # plot(local.BN)
# #
# # # Rows are src nodes and cols are tgt nodes
# # local.BN.adj.matrix <- bnstruct::dag(local.BN)
# # local.BN.adj.matrix <- matrix(local.BN.adj.matrix, nrow = length(local.BN@variables),
# # ncol = length(local.BN@variables),
# # dimnames = c(list(local.BN@variables), list(local.BN@variables)))
# #
# # local.BN.adj.matrix[, centralNodeName]
# # #---------------------------------
# # # End: Local BN struct learning
# # #---------------------------------
#
#
# #---------------------------------
# # Begin: Local Unrolled DBN struct learning
# #---------------------------------
#
# sampleIdx <- 1
#
# local.DBN.input.data <- base::matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
# ncol = base::dim(local.input.data.3D)[2])
#
# for (time.pt in 2:base::dim(local.input.data.3D)[1])
# {
# data.to.combine <- base::matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# ncol = base::dim(local.input.data.3D)[2])
# local.DBN.input.data <- base::cbind(local.DBN.input.data, data.to.combine)
# }
#
# local.DBN.input.data.var.names <- base::c()
# for (time.pt in 1:base::dim(local.input.data.3D)[1]) {
# for (var.name in base::dimnames(local.input.data.3D)[2]) {
# local.DBN.input.data.var.names <- base::c(local.DBN.input.data.var.names,
# base::paste(var.name, as.character(time.pt), sep = "_t"))
# }
# }
# base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
#
# if (base::dim(local.input.data.3D)[3] > 1) # If there are multiple samples per time pt
# {
# for (sampleIdx in 2:base::dim(local.input.data.3D)[3])
# {
# sample.to.combine <- base::matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
# ncol = base::dim(local.input.data.3D)[2])
#
# for (time.pt in 2:base::dim(local.input.data.3D)[1])
# {
# data.to.combine <- base::matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# ncol = base::dim(local.input.data.3D)[2])
# sample.to.combine <- base::cbind(sample.to.combine, data.to.combine)
# }
#
# # local.DBN.input.data.var.names <- c()
# # for (time.pt in 1:dim(local.input.data.3D)[1]) {
# # for (var.name in dimnames(local.input.data.3D)[2]) {
# # local.DBN.input.data.var.names <- c(local.DBN.input.data.var.names,
# # paste(var.name, as.character(time.pt), sep = "_t"))
# # }
# # }
# # colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
#
# local.DBN.input.data <- base::rbind(local.DBN.input.data, sample.to.combine)
# }
# }
#
# base::rownames(local.DBN.input.data) <- base::as.vector(base::unlist(base::dimnames(local.input.data.3D)[3]))
#
# # local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
# # discreteness = rep(TRUE, ncol(local.DBN.input.data)),
# # variables = colnames(local.DBN.input.data),
# # node.sizes = rep(num.discr.levels, ncol(local.DBN.input.data)),
# # starts.from = 0,
# # num.time.steps = nrow(local.input.data.3D))
#
#
# # Begin: Uncomment this section after testing
# local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
# discreteness = base::rep(TRUE, base::ncol(local.DBN.input.data)),
# variables = base::colnames(local.DBN.input.data),
# node.sizes = base::rep(num.discr.levels, base::ncol(local.DBN.input.data)),
# num.time.steps = base::nrow(local.input.data.3D))
#
# # algo = "mmhc", scoring.func = "BDeu", no layering
# # local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# # num.time.steps =
# # bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
#
# # algo = "mmhc", scoring.func = "BDeu", with layering.
# # A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- base::c()
# for (time.pt in 1:num.timepts)
# {
# layers <- base::c(layers, base::rep(time.pt, base::dim(local.input.data.3D)[2]))
# }
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# initial.network = 'random.chain', seed = 12345,
# algo = 'mmhc',
# scoring.func = 'BDeu',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
#
#
# # algo = "sm", scoring.func = "BDeu", no layering
# # local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# # algo = 'sm',
# # num.time.steps =
# # bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
#
#
# # # The following four lines must be executed at the same time
# # save.plot.to.filename = paste(paste('LocalUnrolledDbn', centralNodeName, sep = '_'), '.jpg', sep = '')
# # jpeg(file = paste('LocalUnrolledDbn_', centralNodeName, '.jpg', sep = ''))
# # plot(local.unrolled.DBN)
# # dev.off()
#
# # Extracting the adjacency matrix of the local DBN
# local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
# local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
# nrow = base::length(local.unrolled.DBN@variables),
# ncol = base::length(local.unrolled.DBN@variables),
# dimnames = base::c(base::list(local.unrolled.DBN@variables), base::list(local.unrolled.DBN@variables)))
# # End: Uncomment this section after testing
#
# # # Begin: This section is for testing
# # local.unrolled.DBN.adj.matrix <- matrix(0,
# # nrow = ncol(local.DBN.input.data),
# # ncol = ncol(local.DBN.input.data),
# # dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# # # End: This section is for testing
#
# # 'fixed = FALSE' represents that the given pattern is a regular expression.
# local.unrolled.DBN.central.node.indices <- base::grep(base::paste('^', centralNodeName, sep = ''),
# base::colnames(local.unrolled.DBN.adj.matrix),
# fixed = FALSE)
#
# # Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
#
# #---------------------------------
# # End: Local Unrolled DBN struct learning
# #---------------------------------
#
# # Return value for each 'foreach' iteration
# local.unrolled.DBN.adj.submatrix
# }
#
# # print('foreach ended')
#
# # Shut down the cluster
# parallel::stopCluster(cl)
#
# # print('after stopCluster')
#
# # Begin: Unrolled DBN struct learning
#
# # Begin: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
# unrolled.DBN.adj.matrix.node.names <- base::c()
# for (time.pt in 1:num.timepts) {
# for (node.name in base::dimnames(input.data.discr.3D)[2]) {
# unrolled.DBN.adj.matrix.node.names <- base::c(unrolled.DBN.adj.matrix.node.names,
# base::paste(node.name, as.character(time.pt), sep = "_t"))
# }
# }
#
# unrolled.DBN.adj.matrix <- base::matrix(0, nrow = base::length(unrolled.DBN.adj.matrix.node.names),
# ncol = base::length(unrolled.DBN.adj.matrix.node.names),
# dimnames = base::c(base::list(unrolled.DBN.adj.matrix.node.names),
# base::list(unrolled.DBN.adj.matrix.node.names)))
#
# # End: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
#
# for (list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list))
# {
# # 'submatrix.to.combine'
# submatrix.to.combine <- local.unrolled.DBN.adj.matrix.list[[list.idx]]
#
# unrolled.DBN.adj.matrix[base::rownames(submatrix.to.combine), base::colnames(submatrix.to.combine)] <- submatrix.to.combine
# }
#
# # End: Unrolled DBN struct learning
#
# # Return the final directed net of the current iteration
# return (unrolled.DBN.adj.matrix)
# }
#' Unrolled DBN structure learning with Markov Order 0 and 1.
#'
#' Unrolled DBN structure learning with Markov Order 0 and 1.
#' Candidate parents: The target node itself and its CLR net neighbours at immediately previous and current time pt.
#'
#' @param input.data.discr.3D Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
#' @param mi.net.adj.matrix Adjacency matrix of the mutual information network. Rownames and colnames should be node names.
#' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized for performing the DBN structure learning.
#' @param num.nodes number of nodes
#' @param num.timepts number of timepoints
#' @param output.dirname output directory to store files
#'
#' @return Unrolled DBN adjacency matrix
#'
#' @keywords internal
#' @noRd
#' @importFrom foreach %do% %:% when
learnDbnStructLayer3dParDeg1 <- function(input.data.discr.3D, mi.net.adj.matrix, num.discr.levels, num.nodes, num.timepts, output.dirname="./OUTPUT") {
if(!base::is.array(input.data.discr.3D))
{
base::stop("Error in learnDbnStructLayer3dParDeg1 input.data.discr.3D is not an array")
}
if(!base::is.matrix(mi.net.adj.matrix))
{
base::stop("Error in learnDbnStructLayer3dParDeg1 mi.net.adj.matrix is not a matrix")
}
num.time.trans <- (num.timepts - 1)
no_cores <- num.nodes
cl <- parallel::makeCluster(no_cores, outfile = base::paste(output.dirname, 'outfile.txt', sep = '/' ))
doParallel::registerDoParallel(cl)
centralNodeIdx <- NULL
time.trans.idx <- NULL
# '.verbose = TRUE' is used for debugging
# 'when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0' means when central node does not have any neighbour in the mutual info net
local.unrolled.DBN.adj.matrix.list <-
foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
# when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %:%
foreach::foreach(time.trans.idx = 1:num.time.trans, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
when(base::sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %do% # %dopar%
{
centralNodeName <- base::rownames(mi.net.adj.matrix)[centralNodeIdx]
# List names of the central node's neighbours in mi.net.adj.matrix
nbghNames <- base::c()
if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
{
n <- base::nrow(mi.net.adj.matrix)
for (nbrIdx in 1:n)
{
if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
{
nbghNames <- base::rownames(mi.net.adj.matrix)[nbrIdx]
break
}
}
base::rm(n)
}
else if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
{
nbghNames <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
}
local.net.node.names <- base::c(centralNodeName, nbghNames)
# Select subdataset corr. to the nodes in 'local.net.node.names' and time points 'c(time.trans.idx, (time.trans.idx + 1))'
local.input.data.3D <- input.data.discr.3D[base::c(time.trans.idx, (time.trans.idx + 1)), local.net.node.names, ]
#---------------------------------
# Begin: Local Unrolled DBN struct learning
#---------------------------------
## Begin: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# Say, the local nodes are {v1, v2, v3} and the time points are {t1, t2}.
# Then 'local.DBN.input.data.var.names' contains {v1_t1, v2_t1, v3_t1, v1_t2, v2_t2, v3_t2}.
# In that case, 'local.DBN.input.data' will contain columns corr. to elements in 'local.DBN.input.data.var.names' and
# rows corr. to different samples.
## Begin: Generate first row of 'local.DBN.input.data'
sampleIdx <- 1
local.DBN.input.data <- base::matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
ncol = base::dim(local.input.data.3D)[2])
data.to.combine <- base::matrix(local.input.data.3D[2, , sampleIdx], nrow = 1,
ncol = base::dim(local.input.data.3D)[2])
local.DBN.input.data <- base::cbind(local.DBN.input.data, data.to.combine)
base::rm(data.to.combine)
# for (time.pt in 2:dim(local.input.data.3D)[1])
# {
# data.to.combine <- matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# ncol = dim(local.input.data.3D)[2])
# local.DBN.input.data <- cbind(local.DBN.input.data, data.to.combine)
# }
local.DBN.input.data.var.names <- base::c()
for (time.pt in time.trans.idx:(time.trans.idx + 1)) {
for (var.name in base::dimnames(local.input.data.3D)[2]) {
local.DBN.input.data.var.names <- base::c(local.DBN.input.data.var.names,
base::paste(var.name, base::as.character(time.pt), sep = "_t"))
}
}
base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
## End: Generate first row of 'local.DBN.input.data'
## Begin: Generate rest of the rows of 'local.DBN.input.data'
if (base::dim(local.input.data.3D)[3] > 1) # If there are multiple samples per time pt
{
for (sampleIdx in 2:base::dim(local.input.data.3D)[3])
{
sample.to.combine <- base::matrix(local.input.data.3D[1, , sampleIdx], nrow = 1,
ncol = base::dim(local.input.data.3D)[2])
data.to.combine <- base::matrix(local.input.data.3D[2, , sampleIdx], nrow = 1,
ncol = base::dim(local.input.data.3D)[2])
sample.to.combine <- base::cbind(sample.to.combine, data.to.combine)
# rm(data.to.combine)
# for (time.pt in 2:dim(local.input.data.3D)[1])
# {
# data.to.combine <- matrix(local.input.data.3D[time.pt, , sampleIdx], nrow = 1,
# ncol = dim(local.input.data.3D)[2])
# sample.to.combine <- cbind(sample.to.combine, data.to.combine)
# }
# local.DBN.input.data.var.names <- c()
# for (time.pt in 1:dim(local.input.data.3D)[1]) {
# for (var.name in dimnames(local.input.data.3D)[2]) {
# local.DBN.input.data.var.names <- c(local.DBN.input.data.var.names,
# paste(var.name, as.character(time.pt), sep = "_t"))
# }
# }
# colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
local.DBN.input.data <- base::rbind(local.DBN.input.data, sample.to.combine)
}
}
base::rownames(local.DBN.input.data) <- base::as.vector(base::unlist(base::dimnames(local.input.data.3D)[3]))
## End: Generate rest of the rows of 'local.DBN.input.data'
## End: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
# discreteness = rep(TRUE, ncol(local.DBN.input.data)),
# variables = colnames(local.DBN.input.data),
# node.sizes = rep(num.discr.levels, ncol(local.DBN.input.data)),
# starts.from = 0,
# num.time.steps = nrow(local.input.data.3D))
# Begin: Uncomment this section after testing
local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
discreteness = base::rep(TRUE, base::ncol(local.DBN.input.data)),
variables = base::colnames(local.DBN.input.data),
node.sizes = base::rep(num.discr.levels, base::ncol(local.DBN.input.data)),
num.time.steps = base::nrow(local.input.data.3D))
# algo = "mmhc", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# algo = "sm", scoring.func = "BIC", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# scoring.func = 'BIC',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# algo = "mmhc", scoring.func = "BDeu", with layering, initial.network = 'random.chain', seed = 12345
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
layers <- base::c(base::rep(1, base::dim(local.input.data.3D)[2]), base::rep(2, base::dim(local.input.data.3D)[2]))
local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
initial.network = 'random.chain',
seed = 12345,
algo = 'mmhc',
scoring.func = 'BDeu',
num.time.steps =
bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
layering = layers)
# algo = "mmhc", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# # The following four lines must be executed at the same time
# save.plot.to.filename = paste(paste('LocalUnrolledDbn', centralNodeName, sep = '_'), '.jpg', sep = '')
# jpeg(file = paste('LocalUnrolledDbn_', centralNodeName, '.jpg', sep = ''))
# plot(local.unrolled.DBN)
# dev.off()
# Extracting the adjacency matrix of the local DBN
local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
nrow = base::length(local.unrolled.DBN@variables),
ncol = base::length(local.unrolled.DBN@variables),
dimnames = base::c(base::list(local.unrolled.DBN@variables), base::list(local.unrolled.DBN@variables)))
# End: Uncomment this section after testing
# # Begin: This section is for testing
# local.unrolled.DBN.adj.matrix <- matrix(0,
# nrow = ncol(local.DBN.input.data),
# ncol = ncol(local.DBN.input.data),
# dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# # End: This section is for testing
# 'fixed = FALSE' represents that the given pattern is a regular expression.
# local.unrolled.DBN.central.node.indices <- grep(paste('^', centralNodeName, sep = ''),
# colnames(local.unrolled.DBN.adj.matrix),
# fixed = FALSE)
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
local.unrolled.DBN.src.node.names <- local.unrolled.DBN@variables[1:base::dim(local.input.data.3D)[2]]
local.unrolled.DBN.tgt.node.name <- base::paste(centralNodeName, base::as.character((time.trans.idx + 1)), sep = "_t")
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
local.unrolled.DBN.adj.submatrix <- base::matrix(
local.unrolled.DBN.adj.matrix[local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name],
nrow = base::length(local.unrolled.DBN.src.node.names),
ncol = 1,
dimnames = base::c(base::list(local.unrolled.DBN.src.node.names), base::list(local.unrolled.DBN.tgt.node.name)))
# print(local.unrolled.DBN.adj.submatrix)
#---------------------------------
# End: Local Unrolled DBN struct learning
#---------------------------------
# Return value for each 'foreach' iteration
local.unrolled.DBN.adj.submatrix
}
base::rm(time.trans.idx)
base::rm(centralNodeIdx)
# Shut down the cluster
parallel::stopCluster(cl)
# Begin: Unrolled DBN struct learning
# Begin: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
unrolled.DBN.adj.matrix.node.names <- base::c()
for (time.pt in 1:num.timepts) {
for (node.name in base::dimnames(input.data.discr.3D)[2]) {
unrolled.DBN.adj.matrix.node.names <- base::c(unrolled.DBN.adj.matrix.node.names,
base::paste(node.name, base::as.character(time.pt), sep = "_t"))
}
}
unrolled.DBN.adj.matrix <- base::matrix(0, nrow = base::length(unrolled.DBN.adj.matrix.node.names),
ncol = base::length(unrolled.DBN.adj.matrix.node.names),
dimnames = base::c(base::list(unrolled.DBN.adj.matrix.node.names),
base::list(unrolled.DBN.adj.matrix.node.names)))
# End: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
# 'local.unrolled.DBN.adj.matrix.list' is a list of lists of matrices.
# The outer list contains 'local.unrolled.DBN.adj.matrix.list' inner lists.
# Each inner list contains 'num.time.trans' matrices.
# length(local.unrolled.DBN.adj.matrix.list) = num.nodes
for (outer.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list))
{
# length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = num.time.trans
for (inner.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]))
{
# 'submatrix.to.combine'
submatrix.to.combine <- local.unrolled.DBN.adj.matrix.list[[outer.list.idx]][[inner.list.idx]]
unrolled.DBN.adj.matrix[base::rownames(submatrix.to.combine), base::colnames(submatrix.to.combine)] <- submatrix.to.combine
}
}
# End: Unrolled DBN struct learning
# Return the final directed net of the current iteration
return (unrolled.DBN.adj.matrix)
}
#' Unrolled DBN structure learning with Markov Order 1
#'
#' Candidate parents: The target node itself and its CLR net neighbours at immediately previous time pt.
#'
#' @param input.data.discr.3D Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
#' @param mi.net.adj.matrix Adjacency matrix of the mutual information network. Rownames and colnames should be node names.
#' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized for performing the DBN structure learning.
#' @param num.nodes number of nodes
#' @param num.timepts number of timepoints
#' @param max.fanin maximum incoming edges in the graph
#' @param output.dirname output directory to store files
#'
#' @return Unrolled DBN adjacency matrix
#'
#' @keywords internal
#' @noRd
learnDbnStructMo1Layer3dParDeg1 <- function(input.data.discr.3D, mi.net.adj.matrix, num.discr.levels, num.nodes, num.timepts, max.fanin, output.dirname) {
if(!base::is.array(input.data.discr.3D))
{
base::stop("Error in learnDbnStructMo1Layer3dParDeg1 input.data.discr.3D is not an array")
}
if(!base::is.matrix(mi.net.adj.matrix))
{
base::stop("Error in learnDbnStructMo1Layer3dParDeg1 mi.net.adj.matrix is not a matrix")
}
num.time.trans <- (num.timepts - 1)
## Start and register a parallel backend for parallel computing
## 10 cores to be used for grni server
# no_cores <- min(10, num.nodes, (parallel::detectCores() - 1))
# cl <- parallel::makeCluster(no_cores, outfile = paste(getwd(), 'asset/outfile.txt', sep = '/' ))
# doParallel::registerDoParallel(cl)
# '.verbose = TRUE' is used for debugging
# 'when(sum(mi.net.adj.matrix[, centralNodeIdx]) > 0' means when central node has at least one neighbour in the mutual info net
# Use %do% amd %dopar% for serial and parallel computing, resp.
centralNodeIdx <- NULL
time.trans.idx <- NULL
local.unrolled.DBN.adj.matrix.list <-
foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
# when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %:%
foreach::foreach(time.trans.idx = 1:num.time.trans, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
when(base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 0) %do%
{
centralNodeName <- base::rownames(mi.net.adj.matrix)[centralNodeIdx]
# List names of the central node's neighbours in mi.net.adj.matrix
nbghNames <- base::c()
if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
{
for (nbrIdx in 1:base::nrow(mi.net.adj.matrix))
{
if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
{
nbghNames <- base::rownames(mi.net.adj.matrix)[nbrIdx]
break
}
}
}
else if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
{
nbghNames <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
}
local.net.node.names <- base::c(centralNodeName, nbghNames)
#---------------------------------
# Begin: Local Unrolled DBN struct learning
#---------------------------------
## Begin: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# Say, the central node is v1 and the local nodes are {v1, v2, v3}, and the time points are {t1, t2}.
# Then 'local.DBN.input.data.var.names' contains {v1_t1, v2_t1, v3_t1, v1_t2}.
# In that case, 'local.DBN.input.data' will contain columns corr. to elements in 'local.DBN.input.data.var.names' and
# rows corr. to different samples.
num.samples <- base::dim(input.data.discr.3D)[3]
local.DBN.input.data <- base::matrix(0, nrow = num.samples, ncol = (base::length(local.net.node.names) + 1))
local.unrolled.DBN.src.node.names <- base::c()
for (var.name in local.net.node.names) {
local.unrolled.DBN.src.node.names <- base::c(local.unrolled.DBN.src.node.names,
base::paste(var.name, base::as.character(time.trans.idx), sep = "_t"))
}
local.unrolled.DBN.tgt.node.name <- base::paste(centralNodeName, base::as.character(time.trans.idx + 1), sep = "_t")
local.DBN.input.data.var.names <- base::c(local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name)
base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
for (node.name in local.net.node.names)
{
col.name <- base::paste(node.name, base::as.character(time.trans.idx), sep = "_t")
local.DBN.input.data[, col.name] <- input.data.discr.3D[time.trans.idx, node.name, ]
}
col.name <- base::paste(centralNodeName, base::as.character(time.trans.idx + 1), sep = "_t")
local.DBN.input.data[, col.name] <- input.data.discr.3D[(time.trans.idx + 1), centralNodeName, ]
base::rm(col.name)
## End: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
discreteness = base::rep(TRUE, base::ncol(local.DBN.input.data)),
variables = base::colnames(local.DBN.input.data),
node.sizes = base::rep(num.discr.levels, base::ncol(local.DBN.input.data)))
# algo = "mmhc", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, (ncol(local.DBN.input.data) - 1)), 2)
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# rm(layers)
## algo = "sm", scoring.func = "BIC", with layering
## There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
## The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
## A node with layer idx j can have parents from layer idx i such that i =< j.
layers <- base::c(base::rep(1, (base::ncol(local.DBN.input.data) - 1)), 2)
local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
algo = 'sm',
scoring.func = 'BIC',
layering = layers)
base::rm(layers)
## algo = "sm", scoring.func = "BIC", with layering, with max.fanin
## There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
## The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
## A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, (ncol(local.DBN.input.data) - 1)), 2)
# local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# scoring.func = 'BIC',
# max.fanin = max.fanin,
# layering = layers)
# rm(layers)
# algo = "mmhc", scoring.func = "BDeu", with layering, initial.network = 'random.chain', seed = 12345
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# initial.network = 'random.chain',
# seed = 12345,
# algo = 'mmhc',
# scoring.func = 'BDeu',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# algo = "mmhc", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# # The following four lines must be executed at the same time
# save.plot.to.filename = paste(paste('LocalUnrolledDbn', centralNodeName, sep = '_'), '.jpg', sep = '')
# jpeg(file = paste('LocalUnrolledDbn_', centralNodeName, '.jpg', sep = ''))
# plot(local.unrolled.DBN)
# dev.off()
# Extracting the adjacency matrix of the local DBN
local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
nrow = base::length(local.unrolled.DBN@variables),
ncol = base::length(local.unrolled.DBN@variables),
dimnames = base::c(base::list(local.unrolled.DBN@variables), base::list(local.unrolled.DBN@variables)))
# This for loop checks whether parents are learnt only for 'local.unrolled.DBN.tgt.node.name'. If
# so, then nothing is printed. Otherwise, prints the column(s) corr. to the undesired tgt node(s).
for (col.idx in 1:(base::ncol(local.unrolled.DBN.adj.matrix) - 1))
{
if (base::sum(local.unrolled.DBN.adj.matrix[, col.idx]) > 0)
{
base::print('Erroneous column')
base::print(local.unrolled.DBN.adj.matrix[, col.idx])
}
}
# End: Uncomment this section after testing
# # Begin: This section is for testing
# local.unrolled.DBN.adj.matrix <- matrix(0,
# nrow = ncol(local.DBN.input.data),
# ncol = ncol(local.DBN.input.data),
# dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# # End: This section is for testing
# 'fixed = FALSE' represents that the given pattern is a regular expression.
# local.unrolled.DBN.central.node.indices <- grep(paste('^', centralNodeName, sep = ''),
# colnames(local.unrolled.DBN.adj.matrix),
# fixed = FALSE)
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
local.unrolled.DBN.adj.submatrix <- base::matrix(
local.unrolled.DBN.adj.matrix[local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name],
nrow = base::length(local.unrolled.DBN.src.node.names),
ncol = 1,
dimnames = base::c(base::list(local.unrolled.DBN.src.node.names), base::list(local.unrolled.DBN.tgt.node.name)))
# print(local.unrolled.DBN.adj.submatrix)
#---------------------------------
# End: Local Unrolled DBN struct learning
#---------------------------------
# Return value for each 'foreach' iteration
local.unrolled.DBN.adj.submatrix
}
base::rm(centralNodeIdx)
base::rm(time.trans.idx)
##print('End of foreach loops')
## Shut down the cluster
# parallel::stopCluster(cl)
# Begin: Unrolled DBN struct learning
# Begin: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
unrolled.DBN.adj.matrix.node.names <- base::c()
for (time.pt in 1:num.timepts) {
for (node.name in base::dimnames(input.data.discr.3D)[2]) {
unrolled.DBN.adj.matrix.node.names <- base::c(unrolled.DBN.adj.matrix.node.names,
base::paste(node.name, base::as.character(time.pt), sep = "_t"))
}
}
unrolled.DBN.adj.matrix <- base::matrix(0, nrow = base::length(unrolled.DBN.adj.matrix.node.names),
ncol = base::length(unrolled.DBN.adj.matrix.node.names),
dimnames = base::c(base::list(unrolled.DBN.adj.matrix.node.names),
base::list(unrolled.DBN.adj.matrix.node.names)))
# End: Initialize 'unrolled.DBN.adj.matrix' as a zero matrix
## Debugging
# writeLines('length(local.unrolled.DBN.adj.matrix.list) = ', length(local.unrolled.DBN.adj.matrix.list), '\n')
# print(length(local.unrolled.DBN.adj.matrix.list))
# 'local.unrolled.DBN.adj.matrix.list' is a list of lists of matrices.
# The outer list contains 'local.unrolled.DBN.adj.matrix.list' inner lists.
# Each inner list contains 'num.time.trans' matrices.
# length(local.unrolled.DBN.adj.matrix.list) = num.nodes
for (outer.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list))
{
## Debugging
# writeLines('outer.list.idx = ', outer.list.idx, '\n')
# print(outer.list.idx)
# writeLines('length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = ',
# length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]), '\n')
# print(length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]))
# print(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]])
# if 'local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]' is not an empty list
if (base::length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) > 0)
{
# length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = num.time.trans
for (inner.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]))
{
## Debugging
# writeLines('inner.list.idx = ', inner.list.idx, '\n')
# print(inner.list.idx)
# 'submatrix.to.combine'
submatrix.to.combine <- local.unrolled.DBN.adj.matrix.list[[outer.list.idx]][[inner.list.idx]]
unrolled.DBN.adj.matrix[base::rownames(submatrix.to.combine), base::colnames(submatrix.to.combine)] <- submatrix.to.combine
}
}
}
# End: Unrolled DBN struct learning
# Return the final directed net of the current iteration
return (unrolled.DBN.adj.matrix)
}
#' Unrolled DBN structure learning with Markov Order 1
#'
#' This function is the newer version of learnDbnStructMo1Layer3dParDeg1().
#' The only difference is in the size of unrolled DBN adjacency matrix. In earlier version,
#' the size is ((V \ times T) \ times (v \ times T)) where V = number of nodes and
#' T = number of time points. But the size is too large when V is very large. For
#' e.g., when V = 4028, the function can not execute successfully even with 32 GB main
#' memory in grni server. In this version, the size is reduced by storing the unrolled DBN
#' in an adjacency list of length (T - 1). The t^{th} element in the list is
#' the predicted network adjacency matrix at the t^{th} time interval. Therefore, each list element
#' is a binary matrix of dimension (V \ times V). Hence, the total size of the unrolled DBN
#' adjacency list is ((T - 1) \ times (V \ times V)).
#' Candidate parents: The target node itself and its CLR net neighbours at immediately previous time pt.
#'
#' @param input.data.discr.3D Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
#' @param mi.net.adj.matrix Adjacency matrix of the mutual information network. Rownames and colnames should be node names.
#' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized for performing the DBN structure learning.
#' @param num.nodes number of nodes
#' @param num.timepts number of timepoints
#' @param max.fanin maximum incoming edges in the graph
#' @param node.names name of the nodes
#' @param clr.algo clr algorithm to use
#'
#' @return Unrolled DBN adjacency matrix
#'
#' @keywords internal
#' @noRd
learnDbnStructMo1Layer3dParDeg1_v2 <- function(input.data.discr.3D, mi.net.adj.matrix,
num.discr.levels, num.nodes, num.timepts, max.fanin,
node.names, clr.algo) {
if(!base::is.array(input.data.discr.3D))
{
base::stop("Error in learnDbnStructMo1Layer3dParDeg1_v2 input.data.discr.3D is not an array")
}
if(!base::is.matrix(mi.net.adj.matrix))
{
base::stop("Error in learnDbnStructMo1Layer3dParDeg1_v2 mi.net.adj.matrix is not a matrix")
}
num.time.trans <- (num.timepts - 1)
## Start and register a parallel backend for parallel computing
## 10 cores to be used for grni server
# no_cores <- min(10, num.nodes, (parallel::detectCores() - 1))
# cl <- parallel::makeCluster(no_cores, outfile = paste(getwd(), 'asset/outfile.txt', sep = '/' ))
# doParallel::registerDoParallel(cl)
# '.verbose = TRUE' is used for debugging
# 'when(sum(mi.net.adj.matrix[, centralNodeIdx]) > 0' means when central node has at least one neighbour in the mutual info net
# Use %do% amd %dopar% for serial and parallel computing, resp.
centralNodeIdx <-NULL
time.trans.idx <-NULL
local.unrolled.DBN.adj.matrix.list <-
foreach::foreach(centralNodeIdx = 1:num.nodes, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
# when(sum(mi.net.adj.matrix[, centralNodeIdx]) != 0) %:%
foreach::foreach(time.trans.idx = 1:num.time.trans, .packages = base::c('foreach', 'bnstruct'), .verbose = TRUE) %:%
when(base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 0) %do%
{
centralNodeName <- base::rownames(mi.net.adj.matrix)[centralNodeIdx]
# List names of the central node's neighbours in mi.net.adj.matrix
nbghNames <-base::c()
if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) == 1) # Just 1 neighbour
{
for (nbrIdx in 1:base::nrow(mi.net.adj.matrix))
{
if (mi.net.adj.matrix[nbrIdx, centralNodeIdx] == 1)
{
nbghNames <- base::rownames(mi.net.adj.matrix)[nbrIdx]
break
}
}
}
else if (base::sum(mi.net.adj.matrix[, centralNodeIdx]) > 1) # Multiple neighbours
{
nbghNames <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, centralNodeIdx] == 1),])
}
local.net.node.names <- base::c()
if (clr.algo == 'CLR2.1') {
local.net.node.names <- nbghNames
} else if ((clr.algo == 'CLR') | (clr.algo == 'CLR2')) {
local.net.node.names <- base::c(centralNodeName, nbghNames)
}
#---------------------------------
# Begin: Local Unrolled DBN struct learning
#---------------------------------
## Begin: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# Say, the central node is v1 and the local nodes are {v1, v2, v3}, and the time points are {t1, t2}.
# Then 'local.DBN.input.data.var.names' contains {v1_t1, v2_t1, v3_t1, v1_t2}.
# In that case, 'local.DBN.input.data' will contain columns corr. to elements in 'local.DBN.input.data.var.names' and
# rows corr. to different samples.
num.samples <- base::dim(input.data.discr.3D)[3]
local.DBN.input.data <- base::matrix(0, nrow = num.samples, ncol = (base::length(local.net.node.names) + 1))
local.unrolled.DBN.src.node.names <- base::c()
for (var.name in local.net.node.names) {
local.unrolled.DBN.src.node.names <- base::c(local.unrolled.DBN.src.node.names,
base::paste(var.name, base::as.character(time.trans.idx), sep = "_t"))
}
base::rm(var.name)
local.unrolled.DBN.tgt.node.name <- base::paste(centralNodeName, base::as.character(time.trans.idx + 1), sep = "_t")
local.DBN.input.data.var.names <- base::c(local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name)
base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
for (node.name in local.net.node.names) {
col.name <- base::paste(node.name, base::as.character(time.trans.idx), sep = "_t")
local.DBN.input.data[, col.name] <- input.data.discr.3D[time.trans.idx, node.name, ]
}
base::rm(node.name)
col.name <- base::paste(centralNodeName, base::as.character(time.trans.idx + 1), sep = "_t")
local.DBN.input.data[, col.name] <- input.data.discr.3D[(time.trans.idx + 1), centralNodeName, ]
base::rm(col.name)
## End: Generate 2D 'local.DBN.input.data' from 'local.input.data.3D'
# save(local.DBN.input.data, file = 'local.DBN.input.data.RData')
local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
discreteness = base::rep(TRUE, ncol(local.DBN.input.data)),
variables = base::colnames(local.DBN.input.data),
node.sizes = base::rep(num.discr.levels,
base::ncol(local.DBN.input.data)))
# algo = "mmhc", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", layering = c()
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset))
# algo = "sm", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, (ncol(local.DBN.input.data) - 1)), 2)
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# rm(layers)
## algo = "sm", scoring.func = "BIC", with layering
## There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
## The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
## A node with layer idx j can have parents from layer idx i such that i =< j.
layers <- base::c(base::rep(1, (base::ncol(local.DBN.input.data) - 1)), 2)
local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
algo = 'sm',
scoring.func = 'BIC',
layering = layers)
base::rm(layers)
## algo = "sm", scoring.func = "BIC", with layering, with max.fanin
## There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
## The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
## A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, (ncol(local.DBN.input.data) - 1)), 2)
# local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
# algo = 'sm',
# scoring.func = 'BIC',
# max.fanin = max.fanin,
# layering = layers)
# rm(layers)
# algo = "mmhc", scoring.func = "BDeu", with layering, initial.network = 'random.chain', seed = 12345
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# initial.network = 'random.chain',
# seed = 12345,
# algo = 'mmhc',
# scoring.func = 'BDeu',
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# algo = "mmhc", scoring.func = "BDeu", with layering
# There are two time pts. represented by c(time.trans.idx, time.trans.idx + 1).
# The nodes belonging to these time pts. are labeled with layer idx 1 and 2, resp.
# A node with layer idx j can have parents from layer idx i such that i =< j.
# layers <- c(rep(1, dim(local.input.data.3D)[2]), rep(2, dim(local.input.data.3D)[2]))
# local.unrolled.DBN <- bnstruct::learn.dynamic.network(local.DBN.input.data.BNDataset,
# num.time.steps =
# bnstruct::num.time.steps(local.DBN.input.data.BNDataset),
# layering = layers)
# # The following four lines must be executed at the same time
# save.plot.to.filename = paste(paste('LocalUnrolledDbn', centralNodeName, sep = '_'), '.jpg', sep = '')
# jpeg(file = paste('LocalUnrolledDbn_', centralNodeName, '.jpg', sep = ''))
# plot(local.unrolled.DBN)
# dev.off()
# Extracting the adjacency matrix of the local DBN
local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
nrow = base::length(local.unrolled.DBN@variables),
ncol = base::length(local.unrolled.DBN@variables),
dimnames = base::c(base::list(local.unrolled.DBN@variables), base::list(local.unrolled.DBN@variables)))
# This for loop checks whether parents are learnt only for 'local.unrolled.DBN.tgt.node.name'. If
# so, then nothing is printed. Otherwise, prints the column(s) corr. to the undesired tgt node(s).
for (col.idx in 1:(base::ncol(local.unrolled.DBN.adj.matrix) - 1))
{
if (base::sum(local.unrolled.DBN.adj.matrix[, col.idx]) > 0)
{
base::print('Erroneous column')
base::print(local.unrolled.DBN.adj.matrix[, col.idx])
}
}
# End: Uncomment this section after testing
# # Begin: This section is for testing
# local.unrolled.DBN.adj.matrix <- matrix(0,
# nrow = ncol(local.DBN.input.data),
# ncol = ncol(local.DBN.input.data),
# dimnames = list(colnames(local.DBN.input.data), colnames(local.DBN.input.data)))
# # End: This section is for testing
# 'fixed = FALSE' represents that the given pattern is a regular expression.
# local.unrolled.DBN.central.node.indices <- grep(paste('^', centralNodeName, sep = ''),
# colnames(local.unrolled.DBN.adj.matrix),
# fixed = FALSE)
# Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
# local.unrolled.DBN.adj.submatrix <- local.unrolled.DBN.adj.matrix[, local.unrolled.DBN.central.node.indices]
## Assuming there are > 1 time points. Otherwise, 'local.unrolled.DBN.adj.submatrix' would become a vector.
local.unrolled.DBN.adj.submatrix <- base::matrix(
local.unrolled.DBN.adj.matrix[local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name],
nrow = base::length(local.unrolled.DBN.src.node.names),
ncol = 1,
dimnames = base::c(base::list(local.unrolled.DBN.src.node.names), base::list(local.unrolled.DBN.tgt.node.name)))
# print(local.unrolled.DBN.adj.submatrix)
#---------------------------------
# End: Local Unrolled DBN struct learning
#---------------------------------
# Return value for each 'foreach' iteration
local.unrolled.DBN.adj.submatrix
}
base::rm(centralNodeIdx)
base::rm(time.trans.idx)
##print('End of foreach loops')
# save(local.unrolled.DBN.adj.matrix.list, file = paste(getwd(), 'asset/local.unrolled.DBN.adj.matrix.list.RData', sep = '/'))
## Shut down the cluster
# parallel::stopCluster(cl)
# Begin: Unrolled DBN struct learning
# Initialize the unrolled DBN adjacency list.
## It is a list of length (T - 1) where T = number of time points. The t^{th} element in the list is
## the predicted network adjacency matrix at the t^{th} time interval. Therefore, each list element
## is a binary matrix of dimension (V \times V) where V = number of nodes. Hence, the total size
## of the unrolled DBN adjacency list is ((T - 1) \times (V \times V)).
## Each adjacency matrix is initialized with a zero matrix of dimension (V \times V). Its rows
## corr. to soruce nodes and the columns corr. to target nodes.
unrolled.DBN.adj.matrix.list <- base::list()
for (list.idx in 1:num.time.trans)
{
unrolled.DBN.adj.matrix.list[[list.idx]] <- base::matrix(0, nrow = num.nodes,
ncol = num.nodes,
dimnames =
base::c(base::list(node.names),
base::list(node.names)))
}
base::rm(list.idx)
## 'local.unrolled.DBN.adj.matrix.list' is a list of lists of matrices.
## The outer list contains length(local.unrolled.DBN.adj.matrix.list) number of inner lists.
## The length is \ge 1 and \le num.nodes. It is < num.nodes when there exists some node
## without any neighbour in the mutual info net.
##
## Each inner list contains 'num.time.trans' matrices.
for (outer.list.idx in 1:base::length(local.unrolled.DBN.adj.matrix.list))
{
## Debugging
# writeLines('outer.list.idx = ', outer.list.idx, '\n')
# print(outer.list.idx)
# writeLines('length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = ',
# length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]), '\n')
# print(length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]))
# print(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]])
## if 'local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]' is not an empty list
if (base::length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) > 0)
{
## length(local.unrolled.DBN.adj.matrix.list[[outer.list.idx]]) = num.time.trans
## for any value of outer.list.idx
for (inner.list.idx in 1:num.time.trans)
{
## Debugging
# writeLines('inner.list.idx = ', inner.list.idx, '\n')
# print(inner.list.idx)
## 'submatrix.to.combine'
submatrix.to.combine <- local.unrolled.DBN.adj.matrix.list[[outer.list.idx]][[inner.list.idx]]
## Begin: remove the timestamps from the row and col names of submatrix.to.combine.
## E.g., a row or col name 'G2_t34' will be converted to just 'G2'.
submatrix.to.combine.rownames <- base::c()
for (row.idx in 1:base::nrow(submatrix.to.combine))
{
old.rowname <- base::rownames(submatrix.to.combine)[row.idx]
## Substitute '_t[0-9]+' pattern with empty string in old.rowname
new.rowname <- base::sub('_t[0-9]+', '', old.rowname)
submatrix.to.combine.rownames <- base::c(submatrix.to.combine.rownames, new.rowname)
}
base::rm(row.idx)
base::rownames(submatrix.to.combine) <- submatrix.to.combine.rownames
base::rm(submatrix.to.combine.rownames)
submatrix.to.combine.colnames <- base::c()
for (col.idx in 1:base::ncol(submatrix.to.combine))
{
old.colname <- base::colnames(submatrix.to.combine)[col.idx]
## Replace '_t[0-9]+' pattern with empty string in old.rowname
new.colname <- base::sub('_t[0-9]+', '', old.colname)
submatrix.to.combine.colnames <-base::c(submatrix.to.combine.colnames, new.colname)
}
base::rm(col.idx)
base::colnames(submatrix.to.combine) <- submatrix.to.combine.colnames
base::rm(submatrix.to.combine.colnames)
## End: remove the timestamps from the row and col names of submatrix.to.combine.
# print('submatrix.to.combine')
# print(submatrix.to.combine)
#
# print('inner.list.idx')
# print(inner.list.idx)
# print('unrolled.DBN.adj.matrix.list[[inner.list.idx]]')
# print(unrolled.DBN.adj.matrix.list[[inner.list.idx]])
unrolled.DBN.adj.matrix.list[[inner.list.idx]][base::rownames(submatrix.to.combine),
base::colnames(submatrix.to.combine)] <- submatrix.to.combine
}
}
}
base::rm(outer.list.idx)
# End: Unrolled DBN struct learning
return (unrolled.DBN.adj.matrix.list)
}
#' Learns DBN structure of Markov order 1 where candidate parents are selected using the CLR3 algo
#'
#' It is a serial algorithmic implementation.
#'
#' @param input.data.discr.3D Dimensions {1 = time points, 2 = variables, 3 = samples under the same time point}.
#' @param mi.net.adj.matrix.list.filename List of filenames (relative/absolute path) of adjacency matrices of the mutual information network.
#' @param num.discr.levels If input data is discretized, then number of discrete levels for each variable. Else if input data is continuous, then number of levels in which data needs to be discretized for performing the DBN structure learning.
#' @param num.nodes number of nodes
#' @param num.timepts number of timepoints
#' @param max.fanin maximum incoming edges in the graph
#' @param node.names name of the nodes
#' @param unrolled.DBN.adj.matrix.list List of unrolled DBN adjacency matrices
#'
#' @return Unrolled DBN adjacency matrix
#'
#' @keywords internal
#' @noRd
LearnDbnStructMo1Clr3Ser <- function(input.data.discr.3D, mi.net.adj.matrix.list.filename,
num.discr.levels, num.nodes, num.timepts, max.fanin,
node.names, unrolled.DBN.adj.matrix.list) {
if(!base::is.array(input.data.discr.3D))
{
base::stop("Error in LearnDbnStructMo1Clr3Ser input.data.discr.3D is not an array")
}
if(!base::is.list(mi.net.adj.matrix.list.filename))
{
base::stop("Error in LearnDbnStructMo1Clr3Ser mi.net.adj.matrix.list.filename is not a list")
}
if(!base::is.list(unrolled.DBN.adj.matrix.list))
{
base::stop("Error in LearnDbnStructMo1Clr3Ser unrolled.DBN.adj.matrix.list is not a list")
}
## Here, each 'time.pt.idx' represents time interval
## ('time.pt.idx', ('time.pt.idx' + 1))
for (time.pt.idx in 1:(num.timepts - 1)) {
## Load the adjacency matrices of time-interval-specific
## CLR nets from the given file
mi.net.adj.matrix.list <- NULL
base::load(mi.net.adj.matrix.list.filename)
## Extract the CLR net specific to the current time interval
mi.net.adj.matrix <- mi.net.adj.matrix.list[[time.pt.idx]]
base::rm(mi.net.adj.matrix.list.filename)
for (tgt.node.idx in 1:num.nodes) {
## If the current target node has no candidate parent,
## then skip to the next target node
if (base::sum(mi.net.adj.matrix[, tgt.node.idx]) == 0) {
next
}
## E.g., 'v1_t2'
## Do not use rownames
tgt.node.name <- base::colnames(mi.net.adj.matrix)[tgt.node.idx]
## E.g., 'v1'
tgt.node.base.name <- node.names[tgt.node.idx]
tgt.node.name.as.src <- base::paste(tgt.node.base.name, base::as.character(time.pt.idx), sep = "_t")
## List names of the target node's neighbours in 'mi.net.adj.matrix'
ngbr.names <- base::c()
if (base::sum(mi.net.adj.matrix[, tgt.node.idx]) == 1) {
## Just one neighbour
for (ngbr.idx in 1:base::nrow(mi.net.adj.matrix)) {
if (mi.net.adj.matrix[ngbr.idx, tgt.node.idx] == 1) {
## Do not use colnames
ngbr.names <- base::rownames(mi.net.adj.matrix)[ngbr.idx]
break
}
}
base::rm(ngbr.idx)
} else if (base::sum(mi.net.adj.matrix[, tgt.node.idx]) > 1) {
## Multiple neighbours
## Do not use colnames
ngbr.names <- base::rownames(mi.net.adj.matrix[which(mi.net.adj.matrix[, tgt.node.idx] == 1),])
}
local.net.node.names <- base::c(tgt.node.name.as.src, ngbr.names)
##------------------------------------------------------------
## Begin: Local Unrolled DBN struct learning
##------------------------------------------------------------
num.samples.per.timept <- base::dim(input.data.discr.3D)[3]
local.DBN.input.data <- base::matrix(0, nrow = num.samples.per.timept,
ncol = (base::length(local.net.node.names) + 1))
local.unrolled.DBN.src.node.names <- base::c()
local.unrolled.DBN.src.node.names <- base::c(ngbr.names, tgt.node.name.as.src)
local.unrolled.DBN.src.node.names <- base::sort(local.unrolled.DBN.src.node.names)
base::rm(tgt.node.name.as.src)
local.unrolled.DBN.tgt.node.name <- tgt.node.name
local.DBN.input.data.var.names <- base::c(local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name)
base::colnames(local.DBN.input.data) <- local.DBN.input.data.var.names
base::rm(local.DBN.input.data.var.names)
for (node.name in local.net.node.names) {
## Substitute '_t[0-9]+' pattern with empty string in 'node.name'.
## If 'node.name' = 'v1_t1', then
## 'node.base.name' = 'v1'.
node.base.name <- base::sub('_t[0-9]+', '', node.name)
local.DBN.input.data[, node.name] <- input.data.discr.3D[time.pt.idx, node.base.name, ]
}
base::rm(node.name)
local.DBN.input.data[, tgt.node.name] <- input.data.discr.3D[(time.pt.idx + 1), tgt.node.base.name, ]
local.DBN.input.data.BNDataset <- bnstruct::BNDataset(local.DBN.input.data,
discreteness = base::rep(TRUE, base::ncol(local.DBN.input.data)),
variables = base::colnames(local.DBN.input.data),
node.sizes = base::rep(num.discr.levels,
base::ncol(local.DBN.input.data)))
layers <- base::c(base::rep(1, (base::ncol(local.DBN.input.data) - 1)), 2)
local.unrolled.DBN <- bnstruct::learn.network(local.DBN.input.data.BNDataset,
algo = 'sm',
scoring.func = 'BIC',
layering = layers)
base::rm(layers)
## Extracting the adjacency matrix of the local DBN
local.unrolled.DBN.adj.matrix <- bnstruct::dag(local.unrolled.DBN)
local.unrolled.DBN.adj.matrix <- base::matrix(local.unrolled.DBN.adj.matrix,
nrow = base::length(local.unrolled.DBN@variables),
ncol = base::length(local.unrolled.DBN@variables),
dimnames = base::c(base::list(local.unrolled.DBN@variables),
base::list(local.unrolled.DBN@variables)))
# This for loop checks whether parents are learnt only for 'local.unrolled.DBN.tgt.node.name'. If
# so, then nothing is printed. Otherwise, prints the column(s) corr. to the undesired tgt node(s).
for (col.idx in 1:(base::ncol(local.unrolled.DBN.adj.matrix) - 1)) {
if (base::sum(local.unrolled.DBN.adj.matrix[, col.idx]) > 0) {
base::print('Erroneous column')
base::print(local.unrolled.DBN.adj.matrix[, col.idx])
}
}
base::rm(col.idx)
## Assuming there are > 1 time points.
## Otherwise, 'submatrix.to.combine' would become a vector.
submatrix.to.combine <- base::matrix(
local.unrolled.DBN.adj.matrix[local.unrolled.DBN.src.node.names, local.unrolled.DBN.tgt.node.name],
nrow = base::length(local.unrolled.DBN.src.node.names),
ncol = 1,
dimnames = base::c(base::list(local.unrolled.DBN.src.node.names), base::list(local.unrolled.DBN.tgt.node.name)))
##------------------------------------------------------------
## End: Local Unrolled DBN struct learning
##------------------------------------------------------------
## Begin: remove the timestamps from the row names of
## 'submatrix.to.combine'.
## E.g., a row name 'G2_t34' will be converted to just 'G2'.
submatrix.to.combine.rownames <- base::c()
for (row.idx in 1:base::nrow(submatrix.to.combine)) {
old.rowname <- base::rownames(submatrix.to.combine)[row.idx]
## Substitute '_t[0-9]+' pattern with empty string in old.rowname
new.rowname <- base::sub('_t[0-9]+', '', old.rowname)
submatrix.to.combine.rownames <- base::c(submatrix.to.combine.rownames, new.rowname)
}
base::rm(row.idx)
base::rownames(submatrix.to.combine) <- submatrix.to.combine.rownames
base::rm(submatrix.to.combine.rownames)
## End: remove the timestamps from the row names of submatrix.to.combine.
unrolled.DBN.adj.matrix.list[[time.pt.idx]][rownames(submatrix.to.combine),
tgt.node.base.name] <- submatrix.to.combine
}
base::rm(tgt.node.idx)
}
base::rm(time.pt.idx)
return(unrolled.DBN.adj.matrix.list)
}
Try the TGS package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.