library(bnlearn)
library(Rgraphviz)
library(gRain)

#Formuła na tworzenie sieci połączeń poprzez wskazanie węzłów i ich powiązań warunkowych
bl.alarm <- model2network('[Burglar][Earthquake][Alarm|Burglar:Earthquake][News|Earthquake][Watson|Alarm]')
plot(bl.alarm)

#Alternatywna forma definiowania grafu poprzez wytworzenie powiązań metodą set.arc
bl.dag = empty.graph(nodes = c("B", "E", "N", "A", "W")) 
bl.dag = set.arc(bl.dag, from = "E", to = "N")
bl.dag = set.arc(bl.dag, from = "E", to = "A") 
bl.dag = set.arc(bl.dag, from = "B", to = "A") 
bl.dag = set.arc(bl.dag, from = "A", to = "W") 
bl.dag = set.arc(bl.dag, from = "E", to = "N") 
plot(bl.dag)

#Alternatywna forma definiowania grafu poprzez połączenie wszystkich węzłów za pomocą macierzy powiązań
survey.dag = empty.graph(nodes = c("Bu", "Er", "Ne", "Al", "Wa")) 
arc.set = matrix(c("Bu", "Al", 
                   "Er", "Al", 
                   "Er", "Ne", 
                   "Al", "Wa"),
                 byrow = TRUE, ncol = 2,
                 dimnames = list(NULL, c("from", "to"))) 
arcs(survey.dag) = arc.set
plot(survey.dag)

#usuwanie krawędzi w grafie
modified.alarm <- drop.arc(bl.alarm, "Earthquake", "News") 
plot(modified.alarm)
#zmiana kierunku krawędzi w grafie
modified.alarm <- reverse.arc(modified.alarm , "Alarm", "Burglar") 
plot(modified.alarm)
#tworzenie nowej krawędzi w grafie
modified.alarm <- set.arc(modified.alarm, "Earthquake", "Burglar") 
plot(modified.alarm)

plot(modified.alarm)
modified.alarm <- set.arc(modified.alarm, "Watson", "Earthquake")
# Acykliczność jest domyślnie wymuszona przez wszystkie funkcje tworzące graf
# W konsekwencji powstanie błąd i 
# arc.operations(x = x, from = from, to = to, op = "set", 
#                         %check.cycles = check.cycles, 


bl.alarm
# Wyświetlenie sąsiadów węzła
nbr(bl.alarm, node = 'Alarm')
# Wyświetlenie rodziców węzła
parents(bl.dag, node = "A")
# Wyświetlenie potomków węzła
children(bl.alarm, 'Alarm')
# Wyświetlenie ,,Markov blanket''
mb(bl.alarm, node = 'Alarm')
root.nodes(bl.alarm)
leaf.nodes(bl.alarm)

# Możemy dostosować formę wizualizacji
plot(bl.alarm, highlight = list(nodes='Alarm'))
# z użyciem bibliotegki Rgraphviz
hlight <- list(nodes = c("Earthquake"), arcs = c("Earthquake", "News"), col = "blue", textCol = "grey")
pp <- graphviz.plot(bl.alarm, highlight = hlight, layout = 'circo')
Rgraphviz::renderGraph(pp)


#Możemy wyświetlić wszystkie łuki lub poprosić o ścieżkę między parą węzłów:
arcs(bl.alarm)

path(bl.alarm, from = "Burglar", to = "Watson")
path(bl.alarm, from = "Watson", to = "Burglar")


yn <- c("yes","no")
B <- array(dimnames = list(Burglar = yn), dim = 2, c(0.30,0.70))
E <- array(dimnames = list(Earthquake = yn), dim = 2, c(0.35,0.65))
A <- array(dimnames = list(Alarm = yn, Earthquake = yn, Burglar = yn), 
           dim = c(2, 2, 2),
           c(0.95,0.05,0.90,0.10,0.60,0.40,0.01,0.99))
W <- array(dimnames = list(Watson = yn, Alarm = yn), 
           dim = c(2, 2), c(0.80,0.20,0.40,0.60))
N <- array(dimnames = list(News = yn, Earthquake = yn), 
           dim = c(2, 2), c(0.60,0.40,0.01,0.99))
cpts <- list(Burglar = B, Earthquake = E, Alarm = A, Watson = W, News = N)
bl.alarm.fit = custom.fit(bl.alarm, cpts)

# Co przechowuje nauczony model Sieci Bayesa
bl.alarm.fit
# Co przechowuje nauczony model Sieci Bayesa dla zmiennej Earthquake
bl.alarm.fit$Earthquake         

#Wykresy prawdopodobieństw dla poszczególnych węzłów
bn.fit.barchart(bl.alarm.fit$Earthquake)
bn.fit.barchart(bl.alarm.fit$News)
bn.fit.barchart(bl.alarm.fit$Alarm)


# Wnioskowanie pełne
gr.alarm <- as.grain(bl.alarm.fit)
# Aby obliczyć wszystkie prawdopodobieństwa w sieci należy wykonać kompilację
gr.alarm <- compile(gr.alarm)
gr.alarm

#Można teraz pytać o prawdopodobieństwa każdego węzła w sieci
querygrain(object=gr.alarm, nodes="Earthquake")
querygrain(object=gr.alarm, nodes="Watson")
querygrain(object=gr.alarm, nodes="Alarm")
querygrain(object=gr.alarm, nodes=c('News','Watson'))

map(q.wnb)

no <- rep("no", 5)
nodes <- c('Burglar', 'Earthquake', 'Alarm', 'Watson', 'News') 
gr.alarm <- setEvidence(object=gr.alarm, nodes=nodes, states=no)
gr.alarm
pEvidence(gr.alarm) # Prawdopodobienstwo zdarzenia to

#Powrot do ustawien (bez dowodow)
gr.alarm <- retractEvidence(gr.alarm, nodes) 
gr.alarm

# Pytanie jakie jest prawdopodobienstwo wlamania, gdy Watson dzwoni
?
    # 1. setEvidence
    # 2. querygrain
gr.alarm.Wats <- setEvidence(object=gr.alarm, nodes = c("News", "Watson"), states = c("yes","yes"))
querygrain(object=gr.alarm.Wats, nodes="Burglar", type='joint')