knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
set.seed(1)
library(flexsynth)

dp <- dp_control(epsilon = 1, delta = 1e-6, mechanism = "gaussian")
dp

n <- 1500
real <- data.frame(
  id     = seq_len(n),
  age    = round(rnorm(n, 62, 11)),
  sex    = factor(sample(c("F", "M"), n, TRUE, prob = c(0.45, 0.55))),
  smoker = sample(c(FALSE, TRUE), n, TRUE, prob = c(0.7, 0.3))
)
real$sbp <- round(0.6 * real$age + ifelse(real$smoker, 8, 0) + rnorm(n, 90, 10))
head(real)

df_txt <- data.frame(id = seq_len(400),
                     age = round(rnorm(400, 60, 8)),
                     site = sample(c("north", "south", "east", "west"), 400, TRUE),
                     stringsAsFactors = FALSE)          # site is bare character
dp_txt <- dp_control(epsilon = 6, delta = 1e-6, mechanism = "gaussian",
                     bounds = list(age = c(18, 100)))
synth(df_txt, ~ id, privacy = dp_txt, seed = 1)$privacy

dp <- dp_control(
  epsilon    = 2,
  delta      = 1e-6,
  mechanism  = "gaussian",
  bounds     = list(age = c(18, 100), sbp = c(60, 240))
)
res <- synth(real, structure = ~ id, privacy = dp, seed = 1)
res

syn <- as.data.frame(res)
head(syn)

res$privacy

dp_auto <- dp_control(epsilon = 2, delta = 1e-6, mechanism = "gaussian")
res_auto <- synth(real, structure = ~ id, privacy = dp_auto, seed = 1)
res_auto$privacy

b <- list(age = c(18, 100), sbp = c(60, 240))
dp_tree  <- dp_control(epsilon = 6, mechanism = "laplace", dependence = "tree",
                       bins = 8, bounds = b)
dp_indep <- dp_control(epsilon = 6, mechanism = "laplace", dependence = "independent",
                       bins = 8, bounds = b)
s_tree  <- as.data.frame(synth(real, ~ id, privacy = dp_tree,  seed = 1))
s_indep <- as.data.frame(synth(real, ~ id, privacy = dp_indep, seed = 1))
c(real  = cor(real$age, real$sbp),
  tree  = cor(s_tree$age,  s_tree$sbp),
  indep = cor(s_indep$age, s_indep$sbp))

dp_eff <- dp_control(epsilon = 6, mechanism = "laplace", dependence = "tree",
                     bins = 8, bounds = b, structure_frac = 0.25)
synth(real, ~ id, privacy = dp_eff, seed = 1)$privacy

dp_aim <- dp_control(epsilon = 6, delta = 1e-6, mechanism = "gaussian",
                     select = "adaptive", treewidth = 2,
                     bounds = list(age = c(18, 100), sbp = c(60, 240)))
synth(real, ~ id, privacy = dp_aim, seed = 1)$privacy

dp_anneal <- dp_control(epsilon = 6, delta = 1e-6, mechanism = "gaussian",
                        select = "adaptive", treewidth = 2, anneal = TRUE,
                        bounds = list(age = c(18, 100), sbp = c(60, 240)))
synth(real, ~ id, privacy = dp_anneal, seed = 1)$privacy

dp_bayes <- dp_control(epsilon = 6, delta = 1e-6, mechanism = "gaussian",
                       dependence = "tree", degree = 2,
                       bounds = list(age = c(18, 100), sbp = c(60, 240)))
synth(real, ~ id, privacy = dp_bayes, seed = 1)$privacy

dp_pgm <- dp_control(epsilon = 6, delta = 1e-6, mechanism = "gaussian",
                     dependence = "tree", estimator = "pgm",
                     bounds = list(age = c(18, 100), sbp = c(60, 240)))
synth(real, ~ id, privacy = dp_pgm, seed = 1)$privacy

dp_aim <- dp_control(epsilon = 8, delta = 1e-6, mechanism = "gaussian",
                     select = "aim", treewidth = 2,
                     bounds = list(age = c(18, 100), sbp = c(60, 240)))
synth(real, ~ id, privacy = dp_aim, seed = 1)$privacy

dp_aim_anneal <- dp_control(epsilon = 8, delta = 1e-6, mechanism = "gaussian",
                            select = "aim", treewidth = 2, anneal = TRUE,
                            bounds = list(age = c(18, 100), sbp = c(60, 240)))
synth(real, ~ id, privacy = dp_aim_anneal, seed = 1)$privacy

# The default (scoring = "model") is shown above; this opts back to the cheaper
# one-way-product reference.
dp_aim_indep <- dp_control(epsilon = 8, delta = 1e-6, mechanism = "gaussian",
                           select = "aim", treewidth = 2, scoring = "independence",
                           bounds = list(age = c(18, 100), sbp = c(60, 240)))
synth(real, ~ id, privacy = dp_aim_indep, seed = 1)$privacy

for (e in c(0.5, 2, 8)) {
  d <- dp_control(epsilon = e, mechanism = "laplace",
                  dependence = "independent",
                  bounds = list(age = c(18, 100), sbp = c(60, 240)))
  s <- as.data.frame(synth(real, ~ id, privacy = d, seed = 1))
  cat(sprintf("epsilon = %-3g  mean(sbp): real %.1f  syn %.1f\n",
              e, mean(real$sbp), mean(s$sbp)))
}

long <- do.call(rbind, lapply(1:400, function(i) {
  nv <- sample(2:4, 1); s <- numeric(nv); s[1] <- rnorm(1, 130, 12)
  for (t in seq_len(nv)[-1]) s[t] <- 0.85 * s[t - 1] + 0.15 * 130 + rnorm(1, 0, 5)
  data.frame(id = i, visit = seq_len(nv), sbp = round(s),
             sex = sample(c("F", "M"), 1))
}))
long$sex <- factor(long$sex)

dp_long <- dp_control(epsilon = 8, mechanism = "laplace",
                      max_rows_per_person = 4, bounds = list(sbp = c(60, 240)))
res_long <- synth(long, structure = ~ id / visit, privacy = dp_long, seed = 1)
res_long$privacy

lag1 <- function(d) {
  d <- d[order(d$id, d$visit), ]
  prev <- ave(d$sbp, d$id, FUN = function(x) c(NA, head(x, -1)))
  ok <- !is.na(prev); cor(prev[ok], d$sbp[ok])
}
syn_long <- as.data.frame(res_long)
c(real = lag1(long), synthetic = lag1(syn_long))

res_base <- synth(long, structure = ~ id / visit,
                  privacy = dp_control(epsilon = 8, mechanism = "laplace",
                                       max_rows_per_person = 4,
                                       bounds = list(sbp = c(60, 240)),
                                       baseline = "sex"),
                  seed = 1)

# every synthetic person now has a single sex across their visits
syn_base <- as.data.frame(res_base)
max(tapply(syn_base$sex, syn_base$id, function(s) length(unique(s))))
res_base$privacy   # one fewer transition histogram than the release above

res_ho <- synth(long, structure = ~ id / visit,
                privacy = dp_control(epsilon = 8, mechanism = "laplace",
                                     dependence = "tree", max_rows_per_person = 4,
                                     bounds = list(sbp = c(60, 240)),
                                     baseline = "sex",
                                     transition_order = 2, transition_cross = 1),
                seed = 1)
res_ho$privacy   # note: transitions order 2 + 1 cross-parent, same (eps, delta)

set.seed(1)
patients <- data.frame(
  id  = 1:300,
  age = round(rnorm(300, 60, 11)),
  sex = factor(sample(c("F", "M"), 300, TRUE)))
adm <- do.call(rbind, lapply(patients$id, function(pid) {
  n <- rpois(1, 1.4); if (n == 0) return(NULL)
  data.frame(id = pid, admission_id = seq_len(n), los = 1L + rpois(n, 4))
}))

dp_link <- dp_control(
  epsilon = 4, mechanism = "laplace",
  max_rows_per_person = list(admissions = 6),        # <= 6 admissions per patient
  domain = "public",
  bounds = list(age = c(18, 100), los = c(0, 60)))

res_link <- synth_linked(
  tables     = list(patients = patients, admissions = adm),
  structures = list(patients = ~ id, admissions = ~ id / admission_id),
  keys       = list(patients = "id", admissions = c("id", "admission_id")),
  privacy    = dp_link, seed = 1)
res_link$privacy

check_linkage(res_link)

dp_cross <- dp_control(
  epsilon = 4, mechanism = "laplace", cross_table = TRUE,
  max_rows_per_person = list(admissions = 6),
  domain = "public",
  bounds = list(age = c(18, 100), los = c(0, 60)))

res_cross <- synth_linked(
  tables     = list(patients = patients, admissions = adm),
  structures = list(patients = ~ id, admissions = ~ id / admission_id),
  keys       = list(patients = "id", admissions = c("id", "admission_id")),
  privacy    = dp_cross, seed = 1)
res_cross$privacy

set.seed(1)
visits <- do.call(rbind, lapply(patients$id, function(pid) {
  k  <- 2L + rpois(1, 1.2)
  st <- character(k); st[1] <- sample(c("stable", "worse"), 1)
  for (i in 2:k) st[i] <- if (runif(1) < 0.85) st[i - 1]
                          else setdiff(c("stable", "worse"), st[i - 1])
  data.frame(id = pid, visit_num = seq_len(k),
             status = factor(st, levels = c("stable", "worse")))
}))

dp_longi <- dp_control(
  epsilon = 6, delta = 1e-6, mechanism = "gaussian",
  max_rows_per_person = c(visits = 6), longitudinal = "visits",
  domain = "public", bounds = list(age = c(18, 100)))

res_longi <- synth_linked(
  tables     = list(patients = patients, visits = visits),
  structures = list(patients = ~ id, visits = ~ id / visit_num),
  keys       = list(patients = "id", visits = c("id", "visit_num")),
  privacy    = dp_longi, seed = 1)
res_longi$privacy

dp_both <- dp_control(
  epsilon = 6, delta = 1e-6, mechanism = "gaussian",
  max_rows_per_person = c(visits = 6), longitudinal = "visits",
  cross_table = TRUE, domain = "public", bounds = list(age = c(18, 100)))

res_both <- synth_linked(
  tables     = list(patients = patients, visits = visits),
  structures = list(patients = ~ id, visits = ~ id / visit_num),
  keys       = list(patients = "id", visits = c("id", "visit_num")),
  privacy    = dp_both, seed = 1)
res_both$privacy

dp_order2 <- dp_control(
  epsilon = 6, delta = 1e-6, mechanism = "gaussian",
  max_rows_per_person = c(visits = 6), longitudinal = "visits",
  transition_order = 2, domain = "public", bounds = list(age = c(18, 100)))

res_order2 <- synth_linked(
  tables     = list(patients = patients, visits = visits),
  structures = list(patients = ~ id, visits = ~ id / visit_num),
  keys       = list(patients = "id", visits = c("id", "visit_num")),
  privacy    = dp_order2, seed = 1)
res_order2$privacy

dp_tp <- dp_control(
  epsilon = 6, delta = 1e-6, mechanism = "gaussian",
  max_rows_per_person = c(visits = 6), longitudinal = "visits",
  cross_table = TRUE, transition_parent = 1,
  domain = "public", bounds = list(age = c(18, 100)))

res_tp <- synth_linked(
  tables     = list(patients = patients, visits = visits),
  structures = list(patients = ~ id, visits = ~ id / visit_num),
  keys       = list(patients = "id", visits = c("id", "visit_num")),
  privacy    = dp_tp, seed = 1)
res_tp$privacy