## Package and options
library("dbnR")
library("data.table")
library("bnlearn")
library("visNetwork")
options(prompt = "R> ", continue = "+  ", width = 70,
  useFancyQuotes = FALSE, datatable.print.class = FALSE)


## First DBN example image (Section 2.1)

dbn_ex <- dbnR::generate_random_network_exp(n_vars = 3, size = 2,
  min_mu = -5, max_mu = 5, min_sd = 0.5, max_sd = 2, min_coef = -1,
  max_coef = 1, seed = 42)
attr(dbn_ex$net, "size") <- 2
dbnR::plot_dynamic_network(dbn_ex$net, reverse = TRUE)

## Folding an example dataset (Section 3)

df <- data.frame(X1 = c(3, 6, 4, 9), X2 = c(-1, -2, -3, -4))
df
dbnR::fold_dt(df, size = 2)
dbnR::fold_dt(df, size = 3)

## Filtered folding example (Section 3)

df <- data.frame(X1 = c(3, 6, 4, 9, 8, 2), X2 = c(-1, -2, -3, -4, -5, -6),
  idx = c(1, 1, 1, 2, 2, 2))
df
dbnR::fold_dt(subset(df, select=-idx), size = 2)
dbnR::filtered_fold_dt(df, size = 2, id_var = 'idx')[]

## Example structures of the visualization tool (Section 3.4)

rm_vars <- c("stator_winding", "i_d", "u_q", "i_q", "u_d", "pm")
tmp_dt <- data.table::as.data.table(dbnR::motor)
tmp_dt[, eval(rm_vars) := NULL] # We remove some variables so that the plots are clearer


bn_net <- bnlearn::mmhc(tmp_dt) # BN structure
dbnR::plot_static_network(bn_net)

dbn_net<- dbnR::learn_dbn_struc(tmp_dt, size = 2)
plot(dbn_net)

## Example of forecasting a TS (Section 4.2)

dt <- dbnR::motor
dt_train <- dt[1:2800]
dt_test <- dt[2801:3000]
size <- 2
net <- dbnR::learn_dbn_struc(dt_train, size, method = "dmmhc")
f_dt_train <- dbnR::fold_dt(dt_train, size)
f_dt_test <- dbnR::fold_dt(dt_test, size)
fit <- dbnR::fit_dbn_params(net, f_dt_train)
res_fore <- suppressWarnings(dbnR::forecast_ts(f_dt_test, fit,
  obj_vars = c("pm_t_0"), ini = 50, len = 20))

## Example of using 'bnlearn' functions with 'dbnR' objects (Section 5)
dbn_ex <- dbnR::generate_random_network_exp(n_vars = 3, size = 2,
  min_mu = -5, max_mu = 5, min_sd = 0.5, max_sd = 2, min_coef = -1,
  max_coef = 1, seed = 42)
names(dbn_ex$net$nodes)
dbn_ex$net$arcs

dbn_ex$net <- bnlearn::drop.arc(dbn_ex$net, "X0_t_1", "X0_t_0")
dbn_ex$net$arcs

dbn_ex$net <- bnlearn::remove.node(dbn_ex$net, "X0_t_0")
names(dbn_ex$net$nodes)

## Full example of using 'dbnR' (Section 6)

dt <- dbnR::motor
summary(dt)

### Learning the DBN structure

dt_train <- dt[1:2800]
dt_test <- dt[2801:3000]
size <- 2
f_dt_train <- dbnR::fold_dt(dt_train, size)
f_dt_test <- dbnR::fold_dt(dt_test, size)
print(names(f_dt_train))

t <- Sys.time()
net_dmmhc <- dbnR::learn_dbn_struc(dt_train, size, method = "dmmhc",
  f_dt = f_dt_train)
Sys.time() - t
set.seed(42)
t <- Sys.time()
net_psoho <- dbnR::learn_dbn_struc(dt_train, size, method = "psoho",
  f_dt = f_dt_train, n_it = 10)
Sys.time() - t
t <- Sys.time()
net_nat <- dbnR::learn_dbn_struc(dt_train, size, method = "natPsoho",
  f_dt = f_dt_train, n_it = 10)
Sys.time() - t

plot(net_dmmhc, subset_nodes = c(bnlearn::mb(net_dmmhc, "pm_t_0"), "pm_t_0"))
plot(net_nat, subset_nodes = c(bnlearn::mb(net_nat, "pm_t_0"), "pm_t_0"))
plot(net_psoho, subset_nodes = c(bnlearn::mb(net_psoho, "pm_t_0"), "pm_t_0"))

### Fitting the parameters

fit_dmmhc <- dbnR::fit_dbn_params(net_dmmhc, f_dt_train)
fit_psoho <- dbnR::fit_dbn_params(net_psoho, f_dt_train)
fit_nat <- dbnR::fit_dbn_params(net_nat, f_dt_train)
fit_dmmhc$pm_t_0

fit_nat$pm_t_1
summary(f_dt_train[, pm_t_1])

ev_vars <- names(f_dt_test)[grepl("t_1$", names(f_dt_test))]
ev <- f_dt_test[1, .SD, .SDcols = ev_vars]
ev

### Inference with the models

res <- dbnR::mvn_inference(attr(fit_dmmhc, "mu"),
  attr(fit_dmmhc, "sigma"), evidence = ev)
res

res <- dbnR::predict_dt(fit_dmmhc, f_dt_test, obj_nodes = "pm_t_0")

res <- dbnR::forecast_ts(f_dt_test, fit_dmmhc, obj_vars = "pm_t_0",
  ini = 40, len = 30, mode = "exact")

### Prediction intervals

plot_pred_int_95 <- function(orig, pred, sigma, col){
  u_bound = pred + sigma_p*1.96
  l_bound = pred - sigma_p*1.96
  max_val = max(c(u_bound, l_bound))
  min_val = min(c(u_bound, l_bound))
  x = seq(length(pred))

  plot(ts(orig), ylim = c(min_val, max_val), ylab = col)
  polygon(c(x, rev(x)), c(u_bound, rev(l_bound)),
    col = adjustcolor("purple",alpha.f=0.1), lty = 0)
  lines(pred, col="red")
  lines(u_bound, col="purple")
  lines(l_bound, col="purple")
}

cols <- names(f_dt_test)[-8]
sigma_p <- mvn_inference(attr(fit, 'mu'), attr(fit, 'sigma'),
  f_dt_test[1, .SD, .SDcols = cols])$sigma_p[1]
plot_pred_int_95(res$orig$pm_t_0, res$pred$pm_t_0, sigma_p, "pm_t_0")

### Forecasting without knowing future values

res <- dbnR::forecast_ts(f_dt_test[40], fit_dmmhc, obj_vars = "pm_t_0",
  ini = 1, len = 5, mode = "exact", plot_res = FALSE, print_res = FALSE)
res$pred$pm_t_0


### Using the DBN as a simulator of the real world

f_dt_i <- f_dt_test[40:70]
summary(f_dt_i$motor_speed_t_0)
f_dt_i[1:5, motor_speed_t_0]
f_dt_i[, motor_speed_t_0 := -1.4]
res <- dbnR::forecast_ts(f_dt_i, fit_dmmhc, obj_vars = "pm_t_0",
  ini = 1, len = 30, mode = "exact", prov_ev = "motor_speed_t_0")

f_dt_i[, motor_speed_t_0 := seq(from = -1.4, to = -1.1, by = 0.3 / 30)]
res <- dbnR::forecast_ts(f_dt_i, fit_dmmhc, obj_vars = "pm_t_0",
  ini = 1, len = 30, mode = "exact", prov_ev = "motor_speed_t_0")

### Comparing results between the DBN structures

mae <- function(orig, pred){
  return(sum(abs(orig - pred)) / length(orig))
}

eval_fit <- function(f_dt_test, fit, obj_var = "pm_t_0", len = 20){
  res <- 0
  for(i in 1:(dim(f_dt_test)[1]-len))
    res_fore <- dbnR::forecast_ts(f_dt_test, fit, obj_vars = obj_var,
      ini = i, len = len, plot_res = FALSE,print_res = FALSE)
    res <- res + mae(res_fore$orig[,get(obj_var)],
      res_fore$pred[,get(obj_var)])

  res <- res / (dim(f_dt_test)[1] - len)

  cat(paste0("The final MAE of the model is: ", res, "\n"))

  return(res)
}

res_dmmhc <- eval_fit(f_dt_test, fit_dmmhc)
res_psoho <- eval_fit(f_dt_test, fit_psoho)
res_nat <- eval_fit(f_dt_test, fit_nat)

sessionInfo()