R/0-bnet-base-functions.R

In gherardovarando/Rbnet: Define and Manage Bayesian Networks

#' New bnet object
#' 
#' Create a bnet object with prescribed variables and optionally arcs selected 
#' via adjacency matrix.
#' 
#' @param variables vector of variables names
#' @param adj.mat  adjacency matrix, not symmetric, acyclicity will be checked 
#' for every possible arcs.
#' @return a bnet object, that is a list with components $variables, 
#' \code{bnet$variables} is also a list with component named
#'         as variables, each element represents a node and is composed of 
#'         three elements: \code{name}, \code{parents} and
#'          \code{prob} (initially set to NULL until fitting of the bnet)
#' @export
#' @examples
#' bnet1<-new_bnet(c("a","b"))
#' summary(bnet1)   
#' plot(bnet1)
new_bnet<-function(variables,adj.mat=NULL){
  bn<-list()
  bn$info$model<-"Bayesian Network"
  bn$variables<-list()
  if (!is.null(variables)){
  for (a in variables){
    bn$variables[[a]]<-list()
    bn$variables[[a]]$name<-a
    bn$variables[[a]]$parents<-c()
  }
  }
  bn$info$fitted<-F
  class(bn)<-"bnet"
  if (!is.null(adj.mat)){
    for (i in 1:(dim(adj.mat)[1])){
      for (j in 1:(dim(adj.mat)[2])){
        if ((adj.mat[i,j]!=0)&(i!=j)){
          bn<-add_arc.bnet(bn,from=variables[i],to=variables[j],check = T)
        }
      }
    }
  }
  return(bn)
}

#' Variables of an object
#' 
#' This is a generic functions that returns the variables/nodes of an R object.
#' @param object an R object
#' @return vector of string
#' @export
#' @examples
#' bnet1<-new_bnet(c("a","b","c","d"))
#' variables(bnet1)
variables<-function(object) UseMethod("variables",object)


#' variables bnet
#' 
#' @param object an object \code{bnet}
#' @return vector of string
#' @export
variables.bnet<-function(object){
  return(attributes(object$variables)$names)
}


#' Convert a bnlearn object to bnet object
#' 
#' @param bn a bn object from package bnlearn
#' @return a bnet object
#' @seealso bnet_to_bn
#' @export
#' @examples
#' \dontrun{data(asia,package= "bnlearn")
#' bn<-bnlearn::hc(asia)
#' bnet<-bn_to_bnet(bn)
#' plot(bnet)
#' summary(bnet)}
bn_to_bnet<-function(bn){
  mopbn<-list()
  mopbn$variables<-list()
  for (a in attributes(bn$nodes)$names){
    mopbn$variables[[a]]<-list()
    mopbn$variables[[a]]$name<-a
    mopbn$variables[[a]]$parents<-bn$arcs[,1][bn$arc[,2]==a]
  }
  class(mopbn)<-"bnet"
  mopbn$info$bnlearn<-T
  return(mopbn)
}

#' Convert a bnet object to bnlearn object
#' 
#' @param bnet a bnet object 
#' @return a bn object from package bnlearn
#' @seealso bn_to_bnet
#' @export
bnet_to_bn<-function(bnet){
  maxp<-max(sapply(bnet$variables,FUN = function(x){length(x$parents)}))
  bn<-list()
  class(bn)<-"bn"
  bn$nodes<-lapply(bnet$variables,FUN = function(x){
    l<-list()
    l$parents<-x$parents
  })
  attributes(bn$nodes)$names<-variables.bnet(bnet)
  bn$arcs<-arcs.bnet(bnet)
  return(bn)
}

#' adjacency matrix for a bnet object
#' 
#' @param object a bnet object
#' @return a matrix, the adjacency matrix relative to \code{object}
#' @export
adj_matrix.bnet<-function(object){
  adjM<-matrix(nrow=length(variables.bnet(object)),ncol=length(variables.bnet(object)),0)
  colnames(adjM)<-variables.bnet(object)
  rownames(adjM)<-colnames(adjM)
  for (a in variables.bnet(object)){
    if (!is.null(object$variables[[a]]$parents)){
      adjM[object$variables[[a]]$parents,a]<-1}
  }
  return(adjM)
}

#' Plot a bnet object
#' 
#' @param x a bnet object
#' @param ... compatibility with plot.default
#' @return a plot of a bnet structure, with targets (red) and predictors (green) 
#' @seealso points.bmop
#' @export
#' @examples
#' bnet<-new_bnet(c("A","B"),matrix(nrow=2,c(0,1,0,0)))
#' plot(bnet)
plot.bnet<-function(x,y=NULL,...){
  oldpar<-par()
  parnames<-names(oldpar)
  parnames<-parnames[!parnames %in% c("cin",
                                      "cra",
                                      "csi",
                                      "cxy",
                                      "din",
                                      "page") ]
  attrs<-list()
  attrs$node$fontsize<-"24"
  attrs$edge$arrowsize<-"0.5"
  attrs$edge$style<-"bold"
  g<-graph::graphAM(adjMat = adj_matrix.bnet(x),edgemode = "directed")
  if (is.null(x$targets)){ 
    Rgraphviz::plot(g,attrs=attrs )
    par(oldpar[parnames])
    return(invisible() )}
  nAttrs<-list()
  if (is.null(x$predictors)){
    predictors<-variables.bnet(object = x)
    predictors<-predictors[!(predictors %in% x$targets)]
  }
  else{
    predictors<-x$predictors
  }
  nAttrs$fillcolor<-c(rep(x = "red",times = length(x$targets)),
                      rep("green",times = length(predictors)))
  names(nAttrs$fillcolor)<-c(x$targets,predictors)
  Rgraphviz::plot(g,nodeAttrs=nAttrs,attrs=attrs,... )
  par(oldpar[parnames])
  return(invisible())
}



#' Arcs of an object
#' 
#' @param object 
#' @return a matrix with two columns, each row represent an arc from 
#' the first element to the second.
#' @export
arcs<-function(object){ UseMethod("arcs",object) }

#' Arcs of an object
#' 
#' @param object 
#' @return a matrix with two columns, each row represent an arc from 
#' the first element to the second.
#' @export
#' @examples
#' bnet<-new_bnet(c("A","B"),matrix(nrow=2,c(0,1,0,0)))
#' plot(bnet)
#' arcs(bnet)
arcs.bnet<-function(object){
  maxp<-max(sapply(object$variables,FUN = function(x){length(x$parents)}))
  arcs<-matrix(ncol=2,nrow=length(variables.bnet(object))*maxp)
  i<-1
  for (v in variables.bnet(object)){
    for (b in object$variables[[v]]$parents){
      arcs[i,1]<-b
      arcs[i,2]<-v
      i<-i+1
    }
  }
  return(arcs[!is.na(arcs[,1]),])
}

# leaf, internal function
leaf<-function(arcs,nodes=NULL){
  if (is.null(nodes)){
    nodes==unique(arcs[,1])
  }
  ix<-!(nodes %in% arcs[,2])
  return((nodes[ix])[1])
}

# isacyclic, internal function
isacyclic.arcs<-function(arcs,nodes){
  if (is.null(arcs)){return(TRUE)}  
  if (is.null(dim(arcs))){return(TRUE)}  
  if (dim(arcs)[1]==0){ return(TRUE)}
  if (is.null(nodes)){ return(TRUE)}
  l<-leaf(arcs,nodes)
  if (is.null(l)){ return(FALSE)}
  if (is.na(l)){return(FALSE)}
  return(isacyclic.arcs(arcs[arcs[,1]!=l,],nodes[nodes!=l]))
}

#' check if an object fullfil the acyclic property
#' 
#' @param object an R object 
#' @return logical
#' @export
is.acyclic<-function(object){UseMethod("is.acyclic",object)}

#' isacyclic bnet
#' 
#' @param object a bnet object
#' @return logical
is.acyclic.bnet<-function(object){
  return(isacyclic.arcs(arcs = arcs.bnet(object),
                        nodes = variables.bnet(object)))
}

#' Markov blanket of a bnet object
#' 
#' @param object a bnet object
#' @param nodes a vector of nodes of \code{object}
#' @return a vector with the markov blanket of the nodes in \code{nodes}
#' @export
#' @examples
#' \dontrun{data(asia,package="bnlearn")
#' bn<-bnlearn::hc(asia)
#' bnet<-bn_to_bnet(bn)
#' plot(bnet)
#' summary(bnet)
#' markovblanket.bnet(bnet,names(asia)[1])}
markovblanket.bnet<-function(object,nodes){
  mbtotal<-c()
  for (node in nodes){
    child<-c()
    for (a in variables.bnet(object)){
      if (is.element(node,object$variables[[a]]$parents)){
        child<-unique(c(child,a))
      }  
    }
    mb<-child
    for (a in child){
      mb<-unique(c(mb,object$variables[[a]]$parents))
    }
    mbtotal<-unique(c(mb,mbtotal,object$variables[[node]]$parents))
  }
  return(unique(c(mbtotal,nodes))) 
}

#' Summary of a bnet object
#' 
#' @param object a bnet object
#' @param ... compatibility arguments
#' @export
#' @examples
#' \dontrun{data(gaussian.test,package="bnlearn")
#' bn<-bnlearn::hc(gaussian.test)
#' bnet<-bn_to_bnet(bn)
#' summary(bnet)}
summary.bnet<-function(object,...){
  summ<-list()
  str<-object$variables[[1]]$name
  if (!(length(object$variables[[1]]$parents)==0)){
    
    str<-paste(str,"|",paste(object$variables[[1]]$parents,collapse=":"),sep="")
  }
  for (a in object$variables[-1]){
    str<-paste(str,"][",a$name,sep="")
    if (!(length(a$parents)==0)){
      str<-paste(str,"|",paste(a$parents,collapse=":"),sep="")
    }
  }
  summ$structure<-paste("[",str,"]",sep="")
  summ$nodes<-length(variables.bnet(object))
  summ$arcs<-dim(arcs.bnet(object))[1]
  summ$averagemb<-mean(sapply(object$variables,FUN = function(x){
    return(length(markovblanket.bnet(object = object,nodes = x$name))-1)
  }))
  summ$averageparents<-mean(sapply(object$variables,FUN = function(x){
    return(length(x$parents))
  }))
  summ$info<-object$info
  summ$targets<-object$targets
  summ$predictors<-object$predictors
  class(summ)<-"summary.bnet"
  return(summ)
}

# summary internal function
print.summary.bnet<-function(x,...){
  cat(x$info$model)
  cat("\n")
  if (!is.null(x$targets)){
    cat("Targets:    ",paste(x$targets,collapse = ", "),"\n")
  }
  if (!is.null(x$predictors)){
    cat("Predictors: ",paste(x$predictors,collapse=", "),"\n")
  }
  cat("\n")
  cat("Model structure: \n")
  cat("  ",x$structure)
  cat("\n")
  cat(paste("  nodes:                       ",x$nodes,sep="\t \t"))
  cat("\n")
  cat(paste("  arcs:                        ",x$arcs,sep="\t \t"))
  cat("\n")
  cat(paste("  average markov blanket size: ",x$averagemb,sep="\t \t"))
  cat("\n")
  cat(paste("  average parents per node:    ",x$averageparents,sep="\t \t"))
  cat("\n")
  cat("\n")
  if (!is.null(x$info$wrapper)){
  if (x$info$wrapper$state){
  cat("Wrapper informations: \n")
  cat("  evaluation function: ",x$info$wrapper$evalfun,"\n")
  cat("  evaluation:          ",x$info$wrapper$eval,"\n")
  }
  }
}

#' Print a bnet object
#' 
#' @param x a bnet object
#' @param ... coptaibility arguments
#' @export
print.bnet<-function(x,...){
  print(summary.bnet(x))
  #   cat("BMoP Bayesian network \n \n")
  #   for (a in x$variables){
  #     cat(a$name)
  #     if (!(length(a$parents)==0)){
  #       cat("\t parents: ")
  #       cat(paste(a$parents,collapse=","))
  #       cat("\n")
  #     }
  #     cat(a$mop$knots)
  #     cat(a$mop$ctrpoints)
  #     cat("\n")
  #   }
}

#' child of selected nodes in a bnet object
#' 
#' @param object a bnet object
#' @param nodes a vector of nodes of bnet
#' @return a vector with the names of the child of \code{nodes}
#' @export
child.bnet<-function(object,nodes){
  child<-c()
  for (node in nodes){
   for (n in object$variables){
     if (node %in% n$parents){
       child<-c(child,n$name)
     }
   }
  }
  return(unique(child))
}

#' check if a bnet is fitted
#' 
#' @param object a bnet object
#' @return logical
#' @export
is.fitted.bnet<-function(object){
  if (is.null(object$info$fitted)){ return(FALSE)}
  return(object$info$fitted)
}


#' Check if an object is a bnet object
#' 
#' @param x an R object
#' @return logical
#' @export
is.bnet <- function(x){
  return(inherits(x = x,what = "bnet"))
}

#' Convert object to bnet 
#' 
#' ..when possible (matrix, bnlearn, graph, formulas)
#' @param x a R object
#' @param ... parameter for compatibility reason
#' @return a bnet object
#' @export
#' @examples
#' mat <- rbind(c(0, 0, 1, 1),
#' c(0, 0, 1, 1),
#' c(0, 0, 0, 1),
#' c(0, 0, 0, 0))
#' g1 <- graph::graphAM(adjMat=mat,edgemode = "directed")
#' formula <- Y ~ X1 + X2 + X3 + X4 
#' bnet1 <- as.bnet(g1)
#' bnet2 <- as.bnet(mat)
#' bnet3 <- as.bnet(3*mat)
#' bnet4 <- as.bnet(formula)
as.bnet<-function(x,...){
  
  if (inherits(x,what= "bnet")){
    return(x)
  }
  
  if (inherits(x,what = "bn")){
    return(bn_to_bnet(bn = x))
  }
  
  if (is.matrix(x)){
    if (is.null(colnames(x))){
      var<-as.character(1:(dim(x)[2]))
    }
    else{
      var<-colnames(x)
    }
    return(new_bnet(variables =var,adj.mat = x ))
  }
  
  if (is.data.frame(x)){
    return(new_bnet(variables = names(x),adj.mat = as.matrix(x)))
  }
  
  if (is.character(x)){
   return(new_bnet(variables = x))  
  }
  
  
  
  if (class(x)=="graphNEL"){
    x<-as(x,"graphAM")}
  
  if (class(x)=="graphAM"){
    return(new_bnet(variables= graph::nodes(x),adj.mat =as(x,"matrix") ) )
  }
  
  if (inherits(x,"formula")){
    return(naive_bayes.bnet(formula=x))
  }
  
  if (is.na(x)){
    warning("NA introduced by coercion")
    return(NA)
  }
  
  return(as.bnet(as.matrix(x)))
}

#' @importFrom stats AIC BIC
NULL

#' AIC of bnet 
#' 
#' Computes de Akaike Information Criteria for a bnet object over 
#' a dataset of observations
#' @param object a bnet object
#' @param data the data to compute the log-likelihood
#' @param ... compatibility arguments
#' @param k penalization coefficient (k=2 standard AIC)
#' @return Numerical, the Akaike Information Criteria computed.
#' @seealso \code{\link{BIC.bnet}}
#' @export
AIC.bnet <- function(object, data, ... , k=2 ){
  aa <- 0
  for (a in variables(object)){
    aa <- aa + AIC(object$variables[[a]]$prob , k = k , 
                   data = as.data.frame(
                     data[, c(a, object$variables[[a]]$parents )]) ) 
  }
  return(aa)
}

#' BIC of bnet 
#' 
#' Computes de Bayesian Information Criteria for a bnet object over 
#' a dataset of observations
#' @param object a bnet object
#' @param data the data to compute the log-likelihood
#' @param ... compatibility arguments
#' @return Numerical, the Bayesian Information Criteria computed.
#' @seealso \code{\link{AIC.bnet}}
BIC.bnet <- function(object, data, ... ){
 return(AIC.bnet(object,data,k=log(dim(data)[1])))
}

#'Dimension of a bnet object
#'
#'@param x a bnet object
#'@return positive number, the dimension of \code{x}
#'@export
dim.bnet<-function(x){
  if (!is.fitted.bnet(x)){
    warning("bnet object is not fitted, impossible to compute dimension")
    return(NULL)
  }
  return(sum(sapply(variables(x),function(v){
    return(dim(x$variables[[v]]$prob))
  })))
}



#' Convert bnet into list of nodes and arcs
#' 
#' @param x a bnet object
#' @return a list with at least component \code{nodes} and \code{arcs}
#' @export
as.COSA.bnet<- function(x){
  L<-list()
  L$nodes <- variables(x)
  L$mins <- sapply(x$variables,function(a){
    return(lower(a$prob))
    })
  L$maxs  <- sapply(x$variables,function(a){
    return(upper(a$prob))
  }) 
  L$arcs  <- arcs(x)
  return(L)
}


#' Return marginal x,y to plot
#' 
#' @param node id of the node to plot the marginal
#' @param object bnet object
#' @param N number of points 
#' @param MIN 
#' @export
marginal.bnet<-function(node,object,N=1000,MIN=0){
  bmop<-object$variables[[node]]$marginal
  if (is.null(bmop)){
  return(list(x=c(0,1),y=c(1,1)))  
  }
  x<-lower(bmop)+(0:(N-1))*(upper(bmop)-lower(bmop))/(N-1)
  y<- evaluate(x,bmop,MIN=MIN)
  return(list(x=x,y=y))
 }