K2_AnyR.Rmd

---
jupyter:
  jupytext:
    formats: ipynb,Rmd
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.4.1
  kernelspec:
    display_name: R
    language: R
    name: ir
---

# Bayesian network structure from database


Bayesian networks are networks that can provide insight into probabilistic dependencies that exist among variables in a dataset. In particular we focus on using a Bayesian belief network as a model of probabilistic dependency.

A Bayesian belief-network $B_S$ is a directed acyclic graph in which nodes represent domain variables and the links represent probabilistic dependencies. Variables can be discrete of continuous, but in this work we will focus on discrete variables. 

A Bayesian belief structure , $B_S$ is augmented by conditional probabilities, $B_P$, to form a Bayesian belief network $B=(B_S, B_P)$. For each node $x_i$ there is a conditional probability function that relates the node to its immediate predecessor, called parents and denoted with $\pi_i$. If a node has no parents a prior probability $P(x_i)$ is specified.

The key feature of belief networks is their explicit representation of the conditional indipendence and dependence among events. In particular the joint probability of any instantiation of all $n$ variables in a belief network can be calculated as follows:
$$
P(X_1, \ldots,X_n)=\prod_{i=1}^n P(X_i|\pi_i)
$$
where $X_i$ represent the instantiation of the variable $x_i$ and the $\pi_i$ are the parents of $x_i$.

## The basic model
Let $D$ be a database of cases, $Z$ be the set of variables represented by $D$ and $B_{S_i}$, $B_{S_j}$ be two belief networks structures containing exactly those variables that are in Z. Then:
$$
\frac{P(B_{S_i}|D)}{P(B_{S_j}|D)}\stackrel{prod-rule}{=}\frac{P(B_{S_i},D)}{P(B_{S_j},D)}
$$
We now make some assumptions for computing $P(B_S,D)$ efficiently:

1. The database variables, which we denote with $Z$, are **discrete**. We thus have:
$$
P(B_S, D)=\int_{B_P}P(D|B_s,B_P)f(B_P|B_S)P(B_S)dB_P
$$
where $B_P$ is a vector which denotes the conditional probability assignments associated with belief network structure $B_S$, and $f$ is the conditional probability density function over $B_P$ given $B_S$.
2. Cases occur **independently**, given a belief network model. The Eq \ref{discr} becomes:
$$
P(B_S, D)=\int_{B_P}\left[\prod_{h=1}^mP(C_h|B_S,B_P)\right]f(B_P|B_S)P(B_S)dB_P
$$
where $m$ is the number of cases in $D$ and $C_h$ is the $h$-th case in $D$.
3. There are **no** cases that have variables with **missing values**.
4. The density function $f(B_P|B_S)$ is **uniform**.

We now denote with $w_{ij}$ the $j$-th of $q_i$ unique instantiations of the values of the variables in $\pi_i$, relative to the ordering of the cases in $D$. We then define $v_{ik}$ the $k$-th of the $r_i$ possible instantiations of the variable $x_i$ in $Z$. We lastly define $N_{ijk}$ as the number of casas in $D$ in which the variable $x_i$ has value $v_{ik}$ and $\pi_i$ is instantiated as $w_{ij}$.
Calling:
$$
N_{ij}=\sum_{k=1}^{r_i}N_{ijk}
$$
we finally get:
$$
P(B_S,D)=P(B_S)\prod_{i=1}^n\prod_{j=1}^{q_i} \frac{(r_i-1)!}{(N_{ij}+r_i-1)!}\prod_{k=1}^{r_i}N_{ijk}!
$$

Once we have $P(B,D)$ it is easy to find the conditional probability. Let $Q$ be the set of all those belief-network structures that contain just the variables in set $Z$, then:
$$
P(B_{S_i}|D)=\frac{P(B_{S_i},D)}{\sum_{B_S\in Q}P(B_S,D)}
$$

We can also look at smaller part of the network. Let $G$ be a belief network structure, such that the variables in G are a subset of those in $Z$. Let $R$ be the set of those structures in $Q$ that contains $G$ as a subgraph. Then:
$$
P(G|D)=\frac{\sum_{B_S\in R}P(B_S,D)}{\sum_{B_S\in Q}P(B_S,D)}
$$

## K2 algorithm
However the size and the computational time of finding the structure of a belief network increase exponentially with $n$. So an exact method is not possible: we will apply a greedy search algorithm known as K2. This algorithm begins by making the assumption that a node has no parents. and then adds incrementally that parent whose addition most increase the probability of the resulting structure. We stop the procedure when the addition of no single parent does not increase further the probability.
We define:
$$
g(i,\pi_i)=\prod_{j=1}^{q_i} \frac{(r_i-1)!}{(N_{ij}+r_i-1)!}\prod_{k=1}^{r_i}N_{ijk}!
$$
We also define the function Pred($x_i$) which returns the set of nodes that precede $x_i$ in the node ordering. We thus have the procedure:

- *Input*: a set of $n$ nodes, an ordering of the nodes, an upper bound $u$ to the number of parents a node may have and a database $D$ containing $m$ cases.
- *Output*: For each node a printout of the parents of the node.

```{r}
#Needed packages installation
#install.packages("bnlearn")
#install.packages("bnstruct")
```

```{r}
library("bnlearn")
library("bnstruct")
```

```{r}
# Pseudocode, not to run
for(i in 1:n){  # Cycle over the nodes
    pi.i <- c() # Empty vector
    p.old <- g(i, pi.i) # Function defined above
    OKtoGO <- TRUE
    
    while(OKtoGO & length(pi.i)<u ){
        # let z be the node in Pred(x_i)-pi.i that maximizes g(i, pi.i U {z} )
        z <- which.max( g(i, c(pi.i, pred(x_i)) ) )
        p.new <- g(i, pi.i)
        
        if(p.new >p.old){ # The new configuration is more probable
            p.old <- p.new
            pi.i <- c(pi.i, z) # New parents adding z
        }
        else{  OKtoGO <- FALSE }
        
    }
    cat('Node:', x.i, ' Parents of this node:', pi.i)
}
```

```{r}
# Sample dataset D
# (x1, x2, x3) are binary variables
x1 <- c(1, 1, 0, 1, 0, 0, 1, 0, 1, 0)
x2 <- c(0, 1, 0, 1, 0, 1, 1, 0, 1, 0)
x3 <- c(0, 1, 1, 1, 0, 1, 1, 0, 1, 0)
df <- data.frame(x1, x2, x3)
```

$$ f(i, \pi_i) = \prod_{j=1}^{q_i} \frac{(r_i - 1)!}{(N_{ij} + r_i - 1)!} \prod_{k=1}^{r_i} N_{ijk}! $$


$$ \log f(i, \pi_i) = \sum_{j=1}^{q_i} \Big[ \log (r_i - 1)! - \log (N_{ij} + r_i - 1)!  + \sum_{k=1}^{r_i} \log N_{ijk}! \Big]$$

```{r}
#Precompute all needed (log)factorials (1 to m + r -1)
n <- ncol(df)
m <- nrow(df)
r <- 2 #number of possible values

factorials <- log(c(1, 1:(m+r-1)))

for (i in 4:length(factorials)) { #first 3 are already correct
    factorials[i] <- factorials[i] + factorials[i-1]
}
```

```{r}
f <- function(index, parents, database) { #Rename the variables
    
    #Compute N_{ijk}
    p <- length(parents)
    m <- nrow(database)
    todec <- r ** (0:max(0,p-1)) 

    rowcol <- array(dim = m)
    
    index.k <- database[index] + 1 #(0, 1) -> (1, 2)

    if (p > 0)
        index.j <- as.matrix(database[parents])%*%todec + 1 # 10 -> 2
    else 
        index.j <- rep(1,m)

    
    indexes <- cbind(index.j, index.k) 
    indexes <- indexes[ order(indexes[,1], indexes[,2]),  ]
    
    d <- abs( diff( as.matrix(indexes) ) )
    d1 <- d[,1] + d[,2]
    
    idx.unique <- indexes[c(0,which(d1>0), length(d1)+1), ]
    N_ijk <- diff( c(0,which(d1>0), length(d1)+1) )
    
    N <- data.frame('Row'=idx.unique[,1], 'Col'=idx.unique[,2], 'Value'=N_ijk)
    print(N)
    
    #print(N)           
    #Compute N_{ij}
    Nj <- aggregate(.~Row,data=N,sum)[,c(3)]

    
    #Compute logN - already row summed
    N$Value <- factorials[N$Value+1]  
    logN <- aggregate(.~Row,data=N,sum)[,c(3)]
    

    
    factor.first  <- factorials[r] #r_i-1+1
    factor.second <- factorials[Nj+r] #+r_i -1 +1 (because of the indices)
    factor.third  <- logN

    result <- sum(factor.first - factor.second + factor.third) #factor.first could be brought outside to gain precision (and multiply by qi)

    return(result)
}

#0 is first, 1 is second
```

```{r}
c(1,2,1) + matrix(data=1, ncol=3, nrow=3) #ok broadcasting
```

```{r}
#Compute N_{ijk}
index <- 3
#compute case with parents = 2, 3 with only one loop over database

p <- 1
m <- nrow(database)
todec <- r ** (0:max(0,p-1)) 

rowcol <- array(dim = m)

index.k <- database[index] + 1 #(0, 1) -> (1, 2)

if (p > 0)
    index.j <- as.matrix(database[parents])%*%todec + 1 # 10 -> 2
else 
    index.j <- rep(1,m)


indexes <- cbind(index.j, index.k) 
indexes <- indexes[ order(indexes[,1], indexes[,2]),  ]

d <- abs( diff( as.matrix(indexes) ) )
d1 <- d[,1] + d[,2]

idx.unique <- indexes[c(0,which(d1>0), length(d1)+1), ]
N_ijk <- diff( c(0,which(d1>0), length(d1)+1) )

N <- data.frame('Row'=idx.unique[,1], 'Col'=idx.unique[,2], 'Value'=N_ijk)
print(N)
```

```{r}
f(3, c(1,2), df)
```

```{r}
# Vectorization of the maximization over the possible parents of the probability

# Creating a list of length 1 from the previous parents 
i.par <-  list( i.parents ) 
# Creating a list of length length(i.parents.candidates)
i.cand <- as.list(i.parents.candidates)
# Appending to each element of i.cand the list i.par, and so all the previous parents
i.possible.parents <- mapply(append, i.cand, i.par, SIMPLIFY = FALSE)
# Apply f over all possible parents and then making it a vector
possible.funs <- unlist( lapply(i.possible.parents, f, index=i, database=df) )
log.p.new <- max( possible.funs) # max value
newparent.index <- which.max( possible.funs) # index of the max value
```

```{r}
pred <- function(i, ordering){
    # Parameters:
    # - i: int
    #    index of the variable of which we are searching the precedents
    # - ordering: vector of size n or list of n1 vectors of total size n
    #   If vector must contain a valid permutation of the integers [1, 2 ... n], representing the "prior" ordering of variables,
    #   such that the variable with index ordering[i] can have as possible parents only variables with index ordering[j] with j < i
    #   (i.e. variables that "appear before it" in the provided ordering).
    #   If list must contains the layers of the networks, representing the "prior" ordering of variables,
    #   such that the variable with index ordering[[i]] can have as possible parents only variables in the previous layers, and so
    #   with index ordering[[j]] with j<i. Notice that n1<=n, and if it is equal we come back to the previous case. 
    
    if (class(ordering)=='list'){ # Layer ordering, each element of the list is a layer
            
            for (h in 1:length(ordering)){
               if (i %in% ordering[[h]] ){
                   i.pos <- h
                   break # Once we find i.pos we don't need to look at the other layers
               } 
            }
            if( i.pos == 1){ #If i is in the first layer, there are no parents to check
                next
            } else { # Notice the unlist. We treat all the preceding layers equally
                i.parents.candidates <- unlist( ordering[1:i.pos-1] )
            }
        } else { # Normal ordering, the input is a vector
            
            #Compute vector of candidate parents for i
            i.pos <- which(ordering == i) 
       
            if (i.pos == 1) { #If i is at the start of the ordering, there are no parents to check
                next 
            } else {
                i.parents.candidates <- ordering[1:i.pos-1]
            }
        }
    
    return(i.parents.candidates)
}
```

```{r}
# K2 algorithm implementation

# Parameters:
# - database : data.frame of size (m, n).     #? Since it is always numeric, we could use a matrix also here!
#   Dataset containing m observations of n variables, with no missing values.
#   All values must be integers between 0 and r-1.
# - u : integer
#   Maximum number of parents for each node.
# - ordering : vector of size n or list of n1 vectors of total size n
#   If vector must contain a valid permutation of the integers [1, 2 ... n], representing the "prior" ordering of variables,
#   such that the variable with index ordering[i] can have as possible parents only variables with index ordering[j] with j < i
#   (i.e. variables that "appear before it" in the provided ordering).
#   If list must contains the layers of the networks, representing the "prior" ordering of variables,
#   such that the variable with index ordering[[i]] can have as possible parents only variables in the previous layers, and so
#   with index ordering[[j]] with j<i. Notice that n1<=n, and if it is equal we come back to the previous case. 

# Output:
# adj : matrix of size (n, n)
# Adjacency matrix of the most probable Directed Acyclic Graph found given the evidence in the database.
# adj[i,j] is 1 if there is a connection from node j to node i (i.e. if j is a parent of i)  #? We could take the transpose,
# since it is most common to use i -> j

#TODO:
# - Optimize! (kinda done)
# - Add support for layered ordering (if ordering is a list use layers, if it is vector use current method) (DONE, TO TEST)
# - ? Support for missing values in the database (or we could use bnstruct? See the method in the paper)
# (From the paper: just use every possible value, but it is exponentially complex in the
# number of missing values. Another possibility is to assign the value U to missing
# data, thus treating a variable with r levels and missing data as a variable with r+1 levels)
# - ? Organize code with classes
# - ? Fancy stuff (organize in a R package, use roxygen2 to add documentation, and testthat for unit tests)

K2 <- function(database, u, ordering) {    
    n <- ncol(database)
    adj <- matrix(data = 0, nrow = n, ncol = n)
    
    for (i in 1:n) {
        i.parents <- c()
        log.p.old <- f(i, i.parents, database) #log-Probability of a structure with i disconnected
        
        OKToProceed <- TRUE
        
        # Sostituirei la scelta delle posizioni chiamando una funzione definita fuori
        # prec(i, ordering) per non appesantire troppo il codice
        if (class(ordering)=='list'){ # Layer ordering, each element of the list is a layer
            
            for (h in 1:length(ordering)){
               if (i %in% ordering[[h]] ){
                   i.pos <- h
                   break
               } 
            }
            if( i.pos == 1){
                next
            } else {
                i.parents.candidates <- unlist( ordering[1:i.pos-1] )
            }
        } else { # Normal ordering, the input is a vector
            
            #Compute vector of candidate parents for i
            i.pos <- which(ordering == i) #Position of i in the ordering. Only variables before this position can be parents of i
        
            if (i.pos == 1) { #If i is at the start of the ordering, there are no parents to check
                next 
            } else {
                i.parents.candidates <- ordering[1:i.pos-1]
            }
        }
        
        while (OKToProceed & length(i.parents) < u) {
            log.p.new <- -Inf
            newparent.index <- NA
            
            #Select the node from the candidate list that maximizes the structure's probability (greedy algorithm)
            for (candidate in i.parents.candidates) {
                log.p.candidate <- f(i, c(i.parents, candidate), database)
                
                if (log.p.candidate > log.p.new) {
                    log.p.new <- log.p.candidate
                    newparent.index <- candidate
                }
            }
        
            #Accept the new parent candidate only if it increases the total structure's probability
            if (log.p.new > log.p.old) {
                log.p.old <- log.p.new
                i.parents <- c(i.parents, newparent.index)
                adj[i, newparent.index] = 1
                
                #Remove the added parent from the candidates list
                i.parents.candidates <- i.parents.candidates[-newparent.index]
            } else { #If probability does not increase, stop adding parents to i
                OKToProceed <- FALSE
            }
        }
    }
    
    return(adj)
}
```

```{r}
#Print the structure
structure <- K2(df, 2, c(1,2,3))

print(structure)

#Should be:
# [0, 0, 0]
# [1, 0, 0]
# [0, 1, 0]

for (i in 1:nrow(structure)) {
    cat("Node ", i, " : [")
    
    for (j in 1:ncol(structure)) {
        if (structure[i,j]) {
            cat(j)
        }
    }
    
    cat("]\n")
}

```

```{r}


net <- graph_from_adjacency_matrix(t(structure), mode="directed")
l <- cbind(1:3, rep(0,3))
l <- norm_coords(l, ymin=-1, ymax=0, xmin=-1, xmax=4)

plot(net, layout=l) #for some reason the plot size is too big...
```

```{r}
library('bnlearn')
library('bnstruct')
```

```{r}
datab <- learning.test
n <- ncol(datab)
m <- nrow(datab)

df <- matrix(data=0, nrow = m, ncol = n)
one <- c('b')#, 'wide', 'high')
two <- c('c')
for( i in 1:m){
    for(j in 1:n){
        if(datab[i, j] %in% one) df[i,j] <- 1
        if(datab[i, j] %in% two) df[i,j] <- 2
        
    }
}
df <- data.frame(df)
m <- nrow(df)
r <- 3 #number of possible values


log.fact <- function(m, r){
    # Return the logarithm of the first m+r-1 factorials
    # m: int, number of cases in the dataset
    # r: int, number of possible values a variable can assume
    fact <- log( c(1, 1:(m+r-1)) )

    for (i in 4:length(fact)) { #first 3 are already correct
        fact[i] <- fact[i] + fact[i-1]
    }
    return(fact)
}


factorials <- log.fact(m, r)
```

```{r}
#Print the structure

structure <- K2(df, 3, c(1, 2, 3, 4, 6, 5))

print(structure)

for (i in 1:nrow(structure)) {
    cat("Node ", i, " : [")
    
    for (j in 1:ncol(structure)) {
        if (structure[i,j]) {
            cat(j)
        }
    }
    cat("]\n")
}
```

```{r}
data.alarm <- data.frame(sapply(alarm, as.numeric)-1)

m <- nrow(data.alarm)
r <- 4

factorials <- log.fact(m, r)
```

```{r}
for(n in 1:length(al.name)){
    cat(n, ':', al.name[n], ' ', sep='')
}
```

```{r}
order <- c(18, 20, 3, 9, 15, 23, 36, 21, 19, 17, 16, 30, 11, 35, 22, 29, 
           2, 24, 0, 1, 27, 25, 12, 34, 14, 33, 31, 10, 32, 13, 26, 28, 6, 5, 7, 8, 4)+1

structure <- K2(data.alarm, 4, order)

print(structure)
al.name <- colnames(data.alarm)
for (i in 1:nrow(structure)) {
    cat("Node ", al.name[i], " : [")
    
    for (j in 1:ncol(structure)) {
        if (structure[i,j]) {
            cat(al.name[j], ' ')
        }
    }
    cat("]\n")
}
```

```{r}
?cat
```

```{r}
al_net <- load('alarm.rda')
```

```{r}
names <- sapply(c(1:37), function(x) bn[[x]]$node )
adj <- matrix( data=0, nrow=length(names), ncol=length(names))

i <- 1
for(n in bn){
    par <- n$parents
    for(p in par){
        idx <- which( names==p )
        print(idx)
        adj[i, idx[1]] <- 1
    }
    i <- i+1
}
```

```{r}

```

```{r}
names
```

```{r}
colnames(data.alarm)
```

```{r}
c('HIST')
```

```{r}
#TODO
#---Coding---#
#Evaluate the performance (timing)


#Bnstruct

#---Writing---#
#Write comments for everything
```