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

## -----------------------------------------------------------------------------
# Input documentation:
# named_values: a named numeric vector
# sigma_old, sigma_new: the previous and target diffusion constants
# t_from_to: scaling of time (options: ms->s, s->ms, or none)
convert_prms <- function(named_values,
                         sigma_old = 4,
                         sigma_new = 1,
                         t_from_to = "ms->s") {
  # Some rough input checks
  stopifnot(is.numeric(named_values), is.character(names(named_values)))
  stopifnot(is.numeric(sigma_old), is.numeric(sigma_new))
  t_from_to <- match.arg(t_from_to, choices = c("ms->s", "s->ms", "none"))

  # Internal conversion function (takes a name and value pair, and transforms it)
  internal <- function(name, value) {
    name <- match.arg(
      name,
      choices = c("muc", "b", "non_dec", "sd_non_dec", "tau", "a", "A", "alpha")
    )

    # 1. scale for the diffusion constant
    if (name %in% c("muc", "b", "A")) {
      value <- value * (sigma_new / sigma_old)
    }


    # 2. scale for the time
    # determine the scaling per parameter (assuming s->ms)
    scale <- 1
    if (name %in% c("non_dec", "sd_non_dec", "tau")) scale <- 1000
    if (name %in% c("b", "A")) scale <- sqrt(1000)
    if (name %in% c("muc")) scale <- sqrt(1000) / 1000

    # adapt, depending on the t_from_to argument
    if (t_from_to == "ms->s") scale <- 1 / scale
    if (t_from_to == "none") scale <- 1

    value <- value * scale
  }

  # Apply the internal function to each element
  converted_values <- mapply(FUN = internal, names(named_values), named_values)

  return(converted_values)
}

## -----------------------------------------------------------------------------
dmc_s <- dmc_dm()
prms_solve(dmc_s) # current parameter settings for sigma = 1 and seconds
quants_s <- calc_stats(dmc_s, "quantiles") # calculate predicted quantiles
head(quants_s) # show quantiles

# now the same with new diffusion constant of 4 and time scale in milliseconds
dmc_ms <- dmc_dm()
prms_solve(dmc_ms)["sigma"] <- 4 # new diffusion constant
prms_solve(dmc_ms)["t_max"] <- 3000 # 3000 ms is new max time
prms_solve(dmc_ms)["dt"] <- 1 # 1 ms steps
coef(dmc_ms) <- convert_prms(
  named_values = coef(dmc_ms), # the previous parameters
  sigma_old = 1, # diffusion constants
  sigma_new = 4,
  t_from_to = "s->ms" # how shall the time be scaled?
)
quants_ms <- calc_stats(dmc_ms, "quantiles") # calculate predicted quantiles
head(quants_ms) # show quantiles

## -----------------------------------------------------------------------------
coef(dmc_s)
coef(dmc_ms)

## ----echo = F-----------------------------------------------------------------
# "Unit test" the function

# TEST 1 -> converting twice should lead to the same result as previously
a_model <- dmc_dm(instr = "a ~!")
convert_1 <- convert_prms(
  coef(a_model),
  sigma_old = 1,
  sigma_new = 4,
  t_from_to = "s->ms"
)

convert_2 <- convert_prms(
  convert_1,
  sigma_old = 4,
  sigma_new = 1,
  t_from_to = "ms->s"
)

stopifnot(convert_2 == coef(a_model))


# TEST 2 -> quantiles from above should be very similar
stopifnot(
  abs(quants_ms$Quant_corr - quants_s$Quant_corr * 1000) <= 0.001
)


# TEST 3 -> expectation based on "hand" formula
DMCfun_def <- c(A = 20, tau = 30, muc = 0.5, b = 75, non_dec = 300, sd_non_dec = 30, a = 2, alpha = 3)
exp <- convert_prms(DMCfun_def, sigma_new = 0.1, sigma_old = 4, t_from_to = "ms->s")
stopifnot(abs(exp["A"] - 0.01581139) < 0.0001)
stopifnot(abs(exp["tau"] - 0.030) < 0.0001)
stopifnot(abs(exp["muc"] - 0.3952847) < 0.0001)
stopifnot(abs(exp["b"] - 0.05929271) < 0.0001)
stopifnot(abs(exp["non_dec"] - 0.300) < 0.0001)
stopifnot(abs(exp["sd_non_dec"] - 0.030) < 0.0001)
stopifnot(abs(exp["a"] - 2) < 0.0001)
stopifnot(abs(exp["alpha"] - 3) < 0.0001)


# TEST 4 -> expectation based on "hand" formula (no time scaling)
exp <- convert_prms(DMCfun_def, sigma_new = 0.1, sigma_old = 4, t_from_to = "none")
stopifnot(abs(exp["A"] - 0.5) < 0.0001)
stopifnot(abs(exp["tau"] - 30) < 0.0001)
stopifnot(abs(exp["muc"] - 0.0125) < 0.0001)
stopifnot(abs(exp["b"] - 1.875) < 0.0001)
stopifnot(abs(exp["non_dec"] - 300) < 0.0001)
stopifnot(abs(exp["sd_non_dec"] - 30) < 0.0001)
stopifnot(abs(exp["a"] - 2) < 0.0001)
stopifnot(abs(exp["alpha"] - 3) < 0.0001)