library("bnlearn")
source("include.R")

# Parameters
network = "alarm" # "child", "insurance", "alarm", "mildew", "hailfinder", "munin1", "pigs", "link"
training.samplesize = 5000 # 50 100 200 500 1500 5000
test.samplesize = 5000

# Network loading
alarm.fitted = read.dsc(paste("networks/", network, ".dsc", sep=""))
alarm = bn.net(alarm.fitted)

# Data generation
set.seed(12345) # deterministic random
training = gen.dataset.from.fitted.bn(alarm.fitted, training.samplesize)
test = gen.dataset.from.fitted.bn(alarm.fitted, test.samplesize)

########################################################
##                        H2PC                        ##
########################################################

# Phase I - superstruct
ph1.time = system.time((
  superstruct = hpc(
    x = training,
    test = "mi-h",
    alpha = 0.05,
    optimized = FALSE,
    strict = FALSE,
    undirected = TRUE)
  ))

# Phase II - DAG
ph2.time = system.time((
  dag = tabu(
    x = training,
    start = NULL,
    blacklist = arcs.to.be.added(superstruct$arcs, names(superstruct$nodes)),
    score = "bde",
    tabu = 100,
    max.tabu = 15,
    iss = 10,
    optimized = FALSE)
  ))

# H2PC results
h2pc.models = list(
  superstruct = superstruct,
  dag = dag)
h2pc.measures = c(
  ph1.measures(ph1.time, superstruct),
  ph2.measures(ph2.time, dag))

########################################################
##                        MMHC                        ##
########################################################

# Phase I - superstruct
ph1.time = system.time((
  superstruct = mmpc(
    x = training,
    test = "mi-h",
    alpha = 0.05,
    optimized = FALSE,
    strict = FALSE,
    undirected = TRUE)
))

# Phase II - DAG
ph2.time = system.time((
  dag = tabu(
    x = training,
    start = NULL,
    blacklist = arcs.to.be.added(superstruct$arcs, names(superstruct$nodes)),
    score = "bde",
    tabu = 100,
    max.tabu = 15,
    iss = 10,
    optimized = FALSE)
))

# MMHC results
mmhc.models = list(
  superstruct = superstruct,
  dag = dag)
mmhc.measures = c(
  ph1.measures(ph1.time, superstruct),
  ph2.measures(ph2.time, dag))

########################################################
##                  Increase factors                  ##
########################################################

h2pc.measures / mmhc.measures
