## ----setup, echo = FALSE, message=FALSE---------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(dRiftDM)
set.seed(1014)

re_run = FALSE # if TRUE, runs heavy computations and saves them to inst.. 
# -> otherwise just load them (for speed)

## ----dmc-dm-------------------------------------------------------------------
ddm <- dmc_dm()

## ----chunk-03-----------------------------------------------------------------
print(ddm)

## ----chunk-04-----------------------------------------------------------------
five_traces <- simulate_traces(object = ddm, k = 5)
five_traces

## ----chunk-05-----------------------------------------------------------------
five_traces <- simulate_traces(object = ddm, k = 5, add_x = TRUE)
five_traces

## ----plot-results-------------------------------------------------------------
plot(five_traces, col = c("green", "red"))

## ----plot-results-2-----------------------------------------------------------
exp_behavior <- simulate_traces(object = ddm, k = 1, sigma = 0)
plot(exp_behavior, col = c("green", "red"))

## ----chunk-08-----------------------------------------------------------------
sum_stats <- calc_stats(object = ddm, type = c("cafs", "quantiles"))
sum_stats

## ----plot-results-3-----------------------------------------------------------
plot(sum_stats)

## ----chunk-10-----------------------------------------------------------------
sum_stats <- calc_stats(object = ddm, type = c("densities"))
sum_stats

## ----plot-results-4-----------------------------------------------------------
plot(sum_stats, conds = "comp", xlim = c(0, 1.0))
plot(sum_stats, conds = "incomp", xlim = c(0, 1.0))

## ----chunk-12-----------------------------------------------------------------
coef(ddm)

## ----chunk-13-----------------------------------------------------------------
coef(ddm)["muc"] <- 5
coef(ddm)

## ----chunk-14-----------------------------------------------------------------
coef(ddm, select_unique = FALSE)

## ----chunk-15-----------------------------------------------------------------
prms_solve(ddm)

## ----chunk-16-----------------------------------------------------------------
prms_solve(ddm)["dt"] <- .005
prms_solve(ddm)

## ----chunk-17-----------------------------------------------------------------
check_discretization(ddm)

## ----chunk-18-----------------------------------------------------------------
solver(ddm)

## ----chunk-19-----------------------------------------------------------------
b_coding(ddm)

## ----chunk-20-----------------------------------------------------------------
copy <- ddm # to not change the original model object
b_coding(copy)$column <- "Response"
b_coding(copy)$u_name_value <- c("left" = -1)
b_coding(copy)$l_name_value <- c("right" = 1)

b_coding(copy)

## ----chunk-21-----------------------------------------------------------------
calc_stats(copy, "quantiles")

## ----chunk-22-----------------------------------------------------------------
data <- dRiftDM::dmc_synth_data # some synthetic data suitable for DMC that ships with dRiftDM
# the Cond column matches with conds(ddm).
# The Error column matches b_coding(ddm)
# the RT column is in seconds ;)
head(data)

obs_data(ddm) <- data

## ----chunk-23-----------------------------------------------------------------
extr_data <- obs_data(ddm)
head(extr_data)

## ----chunk-24-----------------------------------------------------------------
summary(ddm)

## ----dmc-dm-2-----------------------------------------------------------------
# get some data (here we use a Simon data set provided by Ulrich et al.)
data <- dRiftDM::ulrich_simon_data
data <- data[data$ID %in% 1:4, ] # just the first four individuals

# get a model (here we again use the pre-built DMC model)
ddm <- dmc_dm()

# the provided data is ready-to-use with DMC
head(data)

# now call the estimation routine
all_fits <- estimate_dm(
  drift_dm_obj = ddm,
  obs_data = data,
  verbose = 1 # prints more information about the underlying optimization run
)

## ----chunk-26-----------------------------------------------------------------
coef(ddm)

## ----fit-model----------------------------------------------------------------
all_fits <- estimate_dm(
  drift_dm_obj = ddm,
  obs_data = data,
  control = list(maxit = 2000), # 2000 iterations; see stats::optim
  verbose = 1 # more information about the underlying optimization run
)

## ----chunk-28-----------------------------------------------------------------
coef(ddm)

## ----fit-model-2--------------------------------------------------------------
l_u <- get_lower_upper(ddm)
print(l_u)

all_fits_nmkb <- estimate_dm(
  drift_dm_obj = ddm,
  obs_data = data,
  optimizer = "nmkb",
  lower = l_u$lower,
  upper = l_u$upper,
  verbose = 1 # more information about the underlying optimization run
)

## ----fit-model-3--------------------------------------------------------------
all_fits_nmkb <- estimate_dm(
  drift_dm_obj = ddm,
  obs_data = data,
  optimizer = "nmkb",
  lower = l_u$lower,
  upper = l_u$upper,
  control = list(maxfeval = 2000), # see dfoptim::nmkb
  verbose = 1 # more information about the underlying optimization run
)

## ----echo = FALSE-------------------------------------------------------------
sys_path <- system.file(package = "dRiftDM")
if (re_run) {
  all_fits_deoptim <- estimate_dm(
    drift_dm_obj = ddm,
    obs_data = data,
    lower = l_u$lower,
    upper = l_u$upper,
    n_cores = 4 # you might want to increase this
  )
  
  # save
  use_directory("inst")
  saveRDS(
    object = all_fits_deoptim,
    file = file.path("inst", "drift_dm_vignette_all_fits_deoptim.rds")
  )
  
} else {
  all_fits_deoptim <- readRDS(
    file = file.path(sys_path, "drift_dm_vignette_all_fits_deoptim.rds")
  )
}


## ----fit-model-4, eval=FALSE--------------------------------------------------
# all_fits_deoptim <- estimate_dm(
#   drift_dm_obj = ddm,
#   obs_data = data,
#   lower = l_u$lower,
#   upper = l_u$upper,
#   n_cores = 1 # you might want to increase this
# )

## ----chunk-32-----------------------------------------------------------------
# Example
class(all_fits_deoptim)
print(all_fits_deoptim)

## ----chunk-33-----------------------------------------------------------------
ddm_est <- estimate_dm(
  drift_dm_obj = ddm,
  obs_data = data[data$ID == 1, ],
  optimizer = "Nelder-Mead",
  control = list(maxit = 2000),
  messaging = FALSE, verbose = 0 # just to reduce the amount of information shown
)
class(ddm_est)
print(ddm_est)

## ----chunk-34-----------------------------------------------------------------
cost_function(ddm) <- "rmse"
print(ddm) # see the end of the following output

## ----fit-model-5--------------------------------------------------------------
fits_agg <- estimate_dm(
  drift_dm_obj = ddm,
  obs_data = data,
  approach = "agg_c",
  optimizer = "nmkb",
  lower = l_u$lower,
  upper = l_u$upper
)

## ----plot-results-5-----------------------------------------------------------
check_fit <- calc_stats(object = all_fits, type = c("cafs", "quantiles"))
plot(check_fit)

## ----plot-results-51----------------------------------------------------------
check_fit <- calc_stats(
  object = all_fits,
  type = c("cafs", "quantiles"),
  level = "group",
  resample = TRUE
)
plot(check_fit)

## ----chunk-37-----------------------------------------------------------------
all_coefs <- coef(all_fits)
print(all_coefs)

## ----plot-results-6-----------------------------------------------------------
hist(all_coefs, bundle_plots = FALSE)

## ----chunk-39-----------------------------------------------------------------
calc_stats(all_fits, type = "fit_stats")

## ----chunk-40-----------------------------------------------------------------
summary(all_fits)

## ----echo = FALSE-------------------------------------------------------------
sys_path <- system.file(package = "dRiftDM")

if (re_run) {
  ddm <- dmc_dm(var_non_dec = FALSE)

  l_u <- get_lower_upper(ddm)
  mcmc_list <- estimate_dm(
    drift_dm_obj = ddm,
    obs_data = data,
    approach = "sep_b",
    lower = l_u$lower,
    upper = l_u$upper,
    verbose = 2,
    n_chains = 40,
    n_cores = 4, # increase to parallelize
    burn_in = 150, # for the demo, usually a bit higher
    samples = 150, # for the demo, usually a bit higher
    seed = 1
  )

  # save the chains (without data to avoid unnec. storage of data)
  use_directory("inst")
  for_save = lapply(mcmc_list, \(x) {
    attr(x, "data_model") <- NULL
    x
  })
  saveRDS(
    object = for_save,
    file = file.path("inst", "drift_dm_vignette_mcmc_list.rds")
  )
} else {
  mcmc_list <- readRDS(
    file = file.path(sys_path, "drift_dm_vignette_mcmc_list.rds")
  )
}

## ----dmc-dm-3, eval = FALSE---------------------------------------------------
# # DMC without variability in the non-decision time;
# # makes it a bit easier to estimate, although it is still
# # difficult, because trial numbers are low in the sample data set :)
# ddm <- dmc_dm(var_non_dec = FALSE)
# 
# # lower and upper boundary for the prior distributions
# l_u <- get_lower_upper(ddm)
# mcmc_list <- estimate_dm(
#   drift_dm_obj = ddm,
#   obs_data = data,
#   approach = "sep_b",
#   lower = l_u$lower,
#   upper = l_u$upper,
#   verbose = 2,
#   n_chains = 40,
#   n_cores = 1, # increase to parallelize
#   burn_in = 150, # for the demo, usually a bit higher
#   samples = 150, # for the demo, usually a bit higher
#   seed = 1
# )

## ----plot-results-7-----------------------------------------------------------
one_mcmc <- mcmc_list[[1]]
one_mcmc
coef(one_mcmc)
plot(one_mcmc, bundle_plots = FALSE)

## ----chunk-43, eval = FALSE---------------------------------------------------
# mcmc_hier <- estimate_dm(
#   drift_dm_obj = ddm,
#   obs_data = data,
#   approach = "hier_b",
#   lower = l_u$lower,
#   upper = l_u$upper,
#   verbose = 2,
#   n_chains = 40,
#   n_cores = 1, # increase for speed-up
#   burn_in = 150, # for the demo, usually way higher
#   samples = 150, # for the demo, usually way higher
#   seed = 1
# )

## ----ratcliff-dm--------------------------------------------------------------
ddm <- ratcliff_dm() # a model for demonstration purpose
sim_1 <- simulate_data(object = ddm, n = 200)
head(sim_1)

## ----chunk-45-----------------------------------------------------------------
sim_2 <- simulate_data(
  object = ddm, n = 200, k = 2,
  lower = c(muc = 1, b = .4, non_dec = .2),
  upper = c(muc = 6, b = .8, non_dec = .4)
)

## ----chunk-46-----------------------------------------------------------------
head(sim_2$synth_data)
head(sim_2$prms)

## ----dmc-dm-4-----------------------------------------------------------------
ddm <- dmc_dm()
traces <- simulate_traces(ddm, k = 2)
# although this object is essentially a list of matrices, the class label ...
class(traces)
print(traces) # ... leads to nicely formatted output; but hides the underlying structure

raw <- unpack_obj(traces) # provides the plain list of matrices
head(t(raw$comp))