# Generated by `rjournal_pdf_article()` using `knitr::purl()`: do not edit by hand
# Please edit simlandr-manuscript.Rmd to modify this file

## ----echo=FALSE---------------------------------------------------------------
knitr::opts_chunk$set(tidy=TRUE, tidy.opts=list(width.cutoff=70), cache=FALSE)
RNGkind("L'Ecuyer-CMRG")


## ----metaphor, fig.show="hold", out.width="100%", fig.cap="Illustration of the ball-and-landscape metaphor commonly used in the field of psychopathology.",fig.alt='Three landscape plots, each with a ball resting in one of two basins. The left basin is labeled maladaptive, the right basin healthy. In the first plot, the maladaptive basin is deeper; in the second, both basins are equally deep; in the third, the healthy basin is deeper.',echo=FALSE,message=FALSE----
knitr::include_graphics("figures/metaphor.png")


## ----diagram, fig.show="hold", out.width="100%", fig.cap="The structure and workflow of simlandr.", fig.alt='A flow chart showing the analysis steps in simlandr, with functions listed under each step.', echo=FALSE,message=FALSE----
knitr::include_graphics("figures/diagram.png")


## ----e1, fig.show="hold", out.width="50%", fig.cap="A graphical illustration of the relationship between the activation levels of the two genes. Solid arrows represent positive relationships (i.e., activation) and dashed arrows represent negative relationships (i.e., inhibition).", fig.alt='Two large circles labeled X1 and X2. Dashed arrows labeled b connect the circles in both directions. Each circle has a solid self-loop labeled a and a dashed self-loop labeled k.', echo=FALSE,message=FALSE----
knitr::include_graphics("figures/diagram-e1.png")


## -----------------------------------------------------------------------------
# Load the package.
library(simlandr)

# Specify the simulation function.
b <- 1
k <- 1
S <- 0.5
n <- 4
lambda <- 0.01

drift_gene <- c(
  rlang::expr(z * x^(!!n) / ((!!S)^(!!n) + x^(!!n)) + (!!b) * (!!S)^(!!n) / ((!!S)^(!!n) + y^(!!n)) - (!!k) * x),
  rlang::expr(z * y^(!!n) / ((!!S)^(!!n) + y^(!!n)) + (!!b) * (!!S)^(!!n) / ((!!S)^(!!n) + x^(!!n)) - (!!k) * y),
  rlang::expr(-(!!lambda) * z)
) |> as.expression()

diffusion_gene <- expression(
  0.2,
  0.2,
  0.2
)


## ----eval=FALSE---------------------------------------------------------------
# # Perform a simulation and save the output.
# set.seed(1614)
# single_output_gene <- sim_SDE(drift = drift_gene, diffusion = diffusion_gene, N = 1e6, M = 10, Dt = 0.1, x0 = c(0, 0, 1), keep_full = FALSE)


## ----echo=FALSE---------------------------------------------------------------
# To save time for building the document, we save the output in a file. The same applies to the following examples.
if (!file.exists("data/single_output_gene.RDS")) {
  set.seed(1614)
  single_output_gene <- sim_SDE(drift = drift_gene, diffusion = diffusion_gene, N = 1e6, M = 10, Dt = 0.1, x0 = c(0, 0, 1), keep_full = FALSE)
  saveRDS(single_output_gene, "data/single_output_gene.RDS")
} else {
  single_output_gene <- readRDS("data/single_output_gene.RDS")
}


## -----------------------------------------------------------------------------
single_output_gene2 <- do.call(rbind, single_output_gene)
single_output_gene2 <- cbind(single_output_gene2[, "X"] - single_output_gene2[, "Y"], single_output_gene2[, "Z"])
colnames(single_output_gene2) <- c("delta_x", "a")


## -----------------------------------------------------------------------------
single_output_gene_mcmc_thin <- as.mcmc.list(lapply(single_output_gene, function(x) x[seq(1, nrow(x), by = 100), ]))


## ----converge-gene, message=FALSE, warning=FALSE, fig.cap="The convergence check result for the simulation of the gene expression model. The variables in different simulation stages did not show distributional differences, indicating that the simulation is long enough to provide a reliable estimation of the steady-state distribution.", fig.alt='Trace and density plots for variables x, y, and z.', out.width="100%"----
plot(single_output_gene_mcmc_thin)


## ----warning=FALSE------------------------------------------------------------
l_single_gene_3d <- 
  make_3d_single(single_output_gene2, 
                 x = "delta_x", y = "a", 
                 lims = c(-1.5, 1.5, 0, 1.5), 
                 Umax = 8)


## ----eval=FALSE---------------------------------------------------------------
# plot(l_single_non_grad_3d)


## ----echo=FALSE,message=FALSE,warning=FALSE,results='hide'--------------------
if(!file.exists("figures/3dstatic_gene.png")) {
  plotly::orca(plot(l_single_gene_3d) |>
  plotly::layout(scene = list(
    aspectmode = "manual", aspectratio = list(x = 1.1, y = 1.1, z = 0.66),
    xaxis = list(range = list(-2, 2)),
    yaxis = list(range = list(0, 1.5)),
    camera = list(
      eye = list(
        x = 0.3, y = -1.5, z = 1.5
      )
    )
  )), file = "figures/3dstatic_gene.png", height = 650, width = 750)
}


## ----3dstaticgene, fig.show="hold", out.width="50%", fig.cap="The 3D landscape (potential value as z-axis) for the gene expression model. The left panel is the plot produced by simlandr; the right panel is the potential landscape obtained analytically by Wang et al. (2008), reproduced with the permission of the authors and in accordance with the journal policy.", fig.alt='Two similar landscape plots, each with three basins.', echo=FALSE,message=FALSE,warning=FALSE----
knitr::include_graphics("figures/3dstatic_gene.png")
knitr::include_graphics("figures/wang2011.png")


## ----results='markup', fig.show='hide', warning=FALSE-------------------------
b_single_gene_3d <- calculate_barrier(l_single_gene_3d,
  start_location_value = c(0, 1.2), end_location_value = c(1, 0.2),
  start_r = 0.3, end_r = 0.3
)

get_barrier_height(b_single_gene_3d)


## ----bsingle3dgene, out.width="50%", fig.cap="The landscape for the gene expression model. The local minima are marked as white dots, the saddle points are marked as red dots, and the MEPs are marked as white lines.", fig.alt='A landscape plot with a white line connecting two white dots, passing through a red dot in the middle.', message=FALSE,warning=FALSE----

plot(l_single_gene_3d, 2) + autolayer(b_single_gene_3d)


## -----------------------------------------------------------------------------
batch_arg_set_gene <- new_arg_set()
batch_arg_set_gene <- batch_arg_set_gene |>
  add_arg_ele(
    arg_name = "parameter", ele_name = "b",
    start = 0.5, end = 1.5, by = 0.5
  ) |>
  add_arg_ele(
    arg_name = "parameter", ele_name = "k",
    start = 0.5, end = 1.5, by = 0.5
  )
batch_grid_gene <- make_arg_grid(batch_arg_set_gene)


## ----eval=FALSE---------------------------------------------------------------
# batch_output_gene <- batch_simulation(batch_grid_gene,
#   sim_fun = function(parameter) {
#     b <- parameter[["b"]]
#     k <- parameter[["k"]]
#     drift_gene <- c(
#       rlang::expr(z * x^(!!n) / ((!!S)^(!!n) + x^(!!n)) + (!!b) * (!!S)^(!!n) / ((!!S)^(!!n) + y^(!!n)) - (!!k) * x),
#       rlang::expr(z * y^(!!n) / ((!!S)^(!!n) + y^(!!n)) + (!!b) * (!!S)^(!!n) / ((!!S)^(!!n) + x^(!!n)) - (!!k) * y),
#       rlang::expr(-(!!lambda) * z)
#     ) |> as.expression()
#     set.seed(1614)
#     single_output_gene <- sim_SDE(drift = drift_gene, diffusion = diffusion_gene, N = 1e6, M = 10, Dt = 0.1, x0 = c(0, 0, 1), keep_full = FALSE)
#     single_output_gene2 <- do.call(rbind, single_output_gene)
#     single_output_gene2 <- cbind(single_output_gene2[, "X"] - single_output_gene2[, "Y"], single_output_gene2[, "Z"])
#     colnames(single_output_gene2) <- c("delta_x", "a")
#     single_output_gene2
#   },
#   bigmemory = TRUE
# )


## ----eval=FALSE---------------------------------------------------------------
# saveRDS(batch_output_gene, "batch_output_gene.RDS")
# batch_output_gene <- readRDS("batch_output_gene.RDS") |> attach_all_matrices()


## ----echo=FALSE---------------------------------------------------------------
if (file.exists("data/batch_output_gene.RDS")) {
  batch_output_gene <- readRDS("data/batch_output_gene.RDS") |> attach_all_matrices()} else {
        batch_output_gene <- batch_simulation(batch_grid_gene,
          sim_fun = function(parameter) {
            b <- parameter[["b"]]
            k <- parameter[["k"]]
            drift_gene <- c(
              rlang::expr(z * x^(!!n) / ((!!S)^(!!n) + x^(!!n)) + (!!b) * (!!S)^(!!n) / ((!!S)^(!!n) + y^(!!n)) - (!!k) * x),
              rlang::expr(z * y^(!!n) / ((!!S)^(!!n) + y^(!!n)) + (!!b) * (!!S)^(!!n) / ((!!S)^(!!n) + x^(!!n)) - (!!k) * y),
              rlang::expr(-(!!lambda) * z)
            ) |> as.expression()
            set.seed(1614)
            single_output_gene <- sim_SDE(drift = drift_gene, diffusion = diffusion_gene, N = 1e6, M = 10, Dt = 0.1, x0 = c(0, 0, 1), keep_full = FALSE)
            single_output_gene2 <- do.call(rbind, single_output_gene)
            single_output_gene2 <- cbind(single_output_gene2[, "X"] - single_output_gene2[, "Y"], single_output_gene2[, "Z"])
            colnames(single_output_gene2) <- c("delta_x", "a")
            single_output_gene2
          },
          bigmemory = TRUE
        )
  saveRDS(batch_output_gene, "data/batch_output_gene.RDS")
  }


## ----warning=FALSE, message=FALSE---------------------------------------------
l_batch_gene_3d <- make_3d_matrix(batch_output_gene,
  x = "delta_x", y = "a", cols = "b", rows = "k",
  lims = c(-5, 5, -0.5, 2), h = 0.005,
  Umax = 8,
  kde_fun = "ks", individual_landscape = TRUE
)


## ----results='markup', fig.show='hide', message=FALSE,warning=FALSE-----------
make_barrier_grid_3d(batch_grid_gene,
  start_location_value = c(0, 1.5), end_location_value = c(1, -0.5),
  start_r = 1, end_r = 1, print_template = TRUE
)

bg_gene <- make_barrier_grid_3d(batch_grid_gene, df = structure(list(start_location_value = list(c(0, 1.5), c(0, 1.5), c(0, 1.5), c(0, 1.5), c(0, 1.5), c(0, 1.5), c(0, 1.5), c(0, 1.5), c(0, 1.5)), start_r = list(c(0.2, 1), c(0.2, 1), c(0.2, 1), c(0.2, 0.5), c(0.2, 0.5), c(0.2, 0.5), c(0.2, 0.3), c(0.2, 0.3), c(0.2, 0.3)), end_location_value = list(
  c(2, 0), c(2, 0), c(2, 0), c(1, 0), c(1, 0), c(1, 0), c(1, 0), c(1, 0), c(1, 0)
), end_r = list(
  c(1, 1), c(1, 1), c(1, 1), c(1, 1), c(1, 1), c(1, 1), c(1, 1), c(1, 1), c(1, 1)
)), row.names = c(NA, -9L), class = c(
  "arg_grid",
  "data.frame"
)))


## -----------------------------------------------------------------------------
b_batch_gene_3d <- calculate_barrier(l_batch_gene_3d,
  bg = bg_gene
)


## ----eval=FALSE---------------------------------------------------------------
# b_batch_gene_3d <- calculate_barrier(l_batch_gene_3d, start_location_value = c(0, 1.5), end_location_value = c(1, 0), start_r = 1, end_r = 1)


## ----bbatch3dgene, out.width="100%", fig.cap="The landscape for the gene expression model of different \\( b \\) and \\( k \\) values. The local minima are marked as white dots, the saddle points are marked as red dots, and the MEPs are marked as white lines.",fig.alt='Nine landscape plots arranged in a 3×3 grid. The x-axis is labeled delta x, the y-axis is labeled a. Rows correspond to k values (0.5, 1, 1.5), and columns to b values (0.5, 1, 1.5).', message=FALSE,warning=FALSE----
plot(l_batch_gene_3d) + autolayer(b_batch_gene_3d)


## ----e2, fig.show="hold", out.width="100%", fig.cap="A graphical illustration of the relationships between several important psychological variables in the panic disorder model. Solid arrows represent positive relationships and dashed arrows represent negative relationships.", fig.alt='Four large circles labeled H, A, PT, and E from left to right. Solid arrows go from A to H, A to PT, PT to A, and PT to E. Dashed arrows go from H to A and from E to PT. The solid arrow from A to PT is labeled AS.', echo=FALSE,message=FALSE----
knitr::include_graphics("figures/diagram-e2.png")


## -----------------------------------------------------------------------------
library(PanicModel)

sim_fun_panic <- function(x0, par) {

  # Change several default parameters
  pars <- pars_default 
  # Increase the noise strength to improve sampling efficiency
  pars$N$lambda_N <- 200
  # Make S constant through the simulation
  pars$TS$r_S_a <- 0
  pars$TS$r_S_e <- 0

  # Specify the initial values of A and PT according to the format requirement by `multi_init_simulation()`, while the other variables use the default initial values.
  initial <- initial_default
  initial$A <- x0[1]
  initial$PT <- x0[2]

  # Specify the value of S according to the format requirement by `batch_simulation()`.
  initial$S <- par$S

  # Extract the simulation output from the result by simPanic(). Only keep the core variables. 
  return(
    as.matrix(
      simPanic(1:5000, initial = initial, parameters = pars)$outmat[, c("A", "PT", "E")]
      )
    )
}


## ----eval=FALSE---------------------------------------------------------------
# future::plan("multisession")
# set.seed(1614, kind = "L'Ecuyer-CMRG")
# single_output_panic <- multi_init_simulation(
#   sim_fun = sim_fun_panic,
#   range_x0 = c(0, 1, 0, 1),
#   R = 4,
#   par = list(S = 0.5)
# )


## ----echo=FALSE---------------------------------------------------------------
if (file.exists("data/single_output_panic.RDS")) {
  single_output_panic <- readRDS("data/single_output_panic.RDS")
} else {
  future::plan("multisession")
  set.seed(1614, kind = "L'Ecuyer-CMRG")
  single_output_panic <- multi_init_simulation(
    sim_fun = sim_fun_panic,
    range_x0 = c(0, 1, 0, 1),
    R = 4,
    par = list(S = 0.5)
  )
  saveRDS(single_output_panic, "data/single_output_panic.RDS")
}


## ----converge-panic, message=FALSE, warning=FALSE, fig.cap="The convergence check result for the simulation of the panic disorder model.", fig.alt='Trace and density plots for variables x, y, and z.', out.width="100%"----
plot(single_output_panic)


## ----3dstaticpanic, fig.show="hold", out.width="50%", fig.cap="The 3D landscape (potential value as color) for the panic disorder model", fig.alt='A landscape plot with two basins.', echo=FALSE,message=FALSE,warning=FALSE----
l_single_panic_3d <- make_3d_single(
  single_output_panic |> window(start = 100), 
  x = "A", y = "PT", h = 0.005, lims = c(-1, 1.5, -0.5, 1.5))
plot(l_single_panic_3d, 2)


## -----------------------------------------------------------------------------
batch_arg_grid_panic <- new_arg_set() |>
  add_arg_ele(arg_name = "par", ele_name = "S", start = 0, end = 1, by = 0.5) |>
  make_arg_grid()


## ----eval=FALSE---------------------------------------------------------------
# future::plan("multisession")
# set.seed(1614, kind = "L'Ecuyer-CMRG")
# batch_output_panic <- batch_simulation(
#   batch_arg_grid_panic,
#   sim_fun = function(par) {
#     multi_init_simulation(
#       sim_fun_panic,
#       range_x0 = c(0, 1, 0, 1),
#       R = 4,
#       par = par
#     ) |> window(start = 100)
#   }
# )


## ----echo=FALSE,message=FALSE-------------------------------------------------
if (file.exists("data/batch_output_panic.RDS")) {
  batch_output_panic <- readRDS("data/batch_output_panic.RDS")
} else {
  future::plan("multisession")
  set.seed(1614, kind = "L'Ecuyer-CMRG")
  batch_output_panic <- batch_simulation(
    batch_arg_grid_panic,
    sim_fun = function(par) {
      multi_init_simulation(
        sim_fun_panic,
        range_x0 = c(0, 1, 0, 1),
        R = 4,
        par = par
      ) |> window(start = 100)
    }
  )
  saveRDS(batch_output_panic, "data/batch_output_panic.RDS")
}


## ----lbatch3dpanic, out.width="100%", fig.cap="The landscape for the panic disorder model of different \\( S \\) values. Two landscapes are shown for different variable combinations, \\( A \\) and \\( PT \\), or \\( A \\) and \\( E \\).", fig.alt='Three landscape plots arranged in a row. The x-axis is labeled A, the y-axis is labeled PT. Columns correspond to S values of 0, 0.5, and 1.', message=FALSE,warning=FALSE----
l_batch_panic_3d <- make_3d_matrix(batch_output_panic, x = "A", y = "PT", cols = "S", h = 0.005, lims = c(-1, 1.5, -0.5, 1.5))
plot(l_batch_panic_3d)

