#**************************************************************************
# PROGRAM TYPE: R function
# PURPOSE: This function is designed to generate a simulation dataset that has the
#          following features
#
#			1. User specified dosing history with the same or different dosing intervals
#        (This scenario applies to a dataset with all subjects receiving the same set of sequential dosing regimens)
#     2. Dosing regimens are read in from a data frame
#        (This scenario applies to a dataset with the subjects receiving only one regimen that can be different between the subjects)
#			2. User specified dosing and observation compartments (multiple observation compartments allowed)
#			3. Adding control points to calculate the AUC invervals between doses (for calculation of CAVG in the NONMEM directly)
#			4. Allow the calculation of CMIN, CMAX and CAVG after each dose
#
# VERSION HISTORY:
# Date        Programmer       Version  Description
# ---------  ---------------- -------  ---------------------------------------
# 17-Nov-2020  Jun Shen         1.0    Final version
#
# MODIFICATION HISTORY:
# DATE       PROGRAMMER            MODIFICATION
# ---------  --------------------  ------------------------------------------
# 5-Apr-2021  Jun Shen             Allow multiple records per ID (time varying) for cov.df
#                                   Allow the dose to be introduced based on the body weight at the time of dosing
#                                   Check dataset for NA or -99
# 4-June-2021 Jun Shen             Allow dose adjustment based on a threshold of body weight change (e.g. >10%)
# 16-Jun-2026 Refactor             Restructured for readability: the repeated event-table
#                                   constructions were factored into builder helpers
#                                   (build.conc / build.dh / build.cmt / build.rate /
#                                   build.cmt.switch[.ss] / merge.covariates), argument
#                                   recycling and the covariate merge were de-duplicated,
#                                   and a latent "condition has length > 1" issue with
#                                   smp.time.ss was fixed. Behaviour is unchanged.
#
# Arguments
# read.dose.info = logical value to indicate if the dosing information may be read in through cov.df
# di = Dose intervals (can be a vector if multiple intervals are used)
# num.dose = the number of doses to be simulated (can be a vector), paired with "di"
# dose = the nominal dose (can be a vector), paired with "di"
# dose.unit = the dose unit (e.g. 'mg' or 'mg/kg') (can be a vector), paired with "di"
# dose.cmt = the dosing compartment number, paired with "di"
# sec.dose.cmt = a second dosing compartment, paried with "di". This is NOT for different doses given through different routes. This is for a dose to be introduced into two compartments simultaneously (e.g. an oral dose to be modeled as zero order and first order absorption)
# obs.cmt = observation compartment number, NOT paired with "di", multiple obs.cmt values meant for observations from multiple compartments at each time point
# smp.time = sampling times (vector), paried with "di"
# smp.time.as.is = TRUE/FALSE, the default behavior is to insert a smp.time vector after each dose. When smp.time.as.is is set to TRUE, the smp.time vector will be added once as is
# inf.time = the infusion time (in hours) for i.v. infusion, paried with "di". If a negative value is provided (e.g. -1 or -2), the value will be passed directly to the RATE parameter
# sec.inf.time = this is used to provide a value to the RATE parameter when sec.dose.cmt argument is used (e.g. -2 to indicate to model the duration of infusion) )
# auc.interval = logical value to indicate if auc interval needs to be calculated, this argument can also be a vector to be paired with "di"
# auc.cmt = the AUC compartment number(s) (could be a single value of a vector)
# cov.df = a data frame with covariate values. One record per ID. Minimally can be just a data frame with one column for ID variable
# id.var = the ID variable (default to "ID") in the cov.df
# bw.var = the body weight variable (default to "BBWT") provide a baseline value for dose calculation
# time.var = The time variable in the cov.df for time varying covariates
# tv.bw.var = The time varying body weight (or any other time varying variable) to be used for calculating the dose based on the body weight at the time of dosing
# bw.change.threshold = a threshold value for the change of body weight above which the dose will be adjusted (e.g. 0.1 means >10% body weight change)
# output.suffix = a string of character used as suffix to be added to the name of the dataset prepared (e.g. 's1')
# export = logical value to indicate if the dataset is exported
# export.path = a path for the dataset to be exported (default to '../../nm/0data/')
# ss = logical value to indicate if a dataset for steady state is needed
# ss.with.addl = Instead of using SS to achieve thearectical steady state, use ADDL to introduce enough doses to achieve equivalent status. this may be useful to compute steady state AUC for non-linear PK
# addl = the additional number of doses to be introduced for steady state (This feature is most likely to be used for nonlinear system when using "SS" item is not appropriate)
# auc.ss.interval = TRUE/FALSE to calculate the auc interval at steady state
# only.ss = logical value to indicate if only the steady state dataset is produced
# smp.time.ss = optional sampling times used for steady state dataset
# num.dose.ss = optional the number of doses to be simulated at steady state (e.g. 2 doses simulated for Q2W)
# skip.dir.check = logical value to indicate if directory check needs to be turned off (due to old system under ValidR)


# ----------------------------------------------------------------------------
# Small utility functions used internally
# ----------------------------------------------------------------------------

# Move character (non-numeric) columns to the end of a data frame
char.col.at.last <- function(data, col.first = NULL) {
  if (!is.data.frame(data)) stop('a data frame is required for this function')
  if (!is.null(col.first)) {
    temp  <- match(col.first, names(data))
    temp2 <- setdiff(1:dim(data)[2], temp)
    data  <- data[c(temp, temp2)]
  }

  char.col <- grep(FALSE, sapply(data, is.numeric))
  num.col  <- grep(TRUE,  sapply(data, is.numeric))

  data <- data[c(num.col, char.col)]
  return(data)
}

# Map values in x from one set of codes to another (NA where unmatched, unless strict=FALSE)
map <- function(x, from, to, strict = TRUE, ...) {
  stopifnot(length(to) == length(from))
  res <- to[match(x, table = from)]
  if (!strict)
    res[!(x %in% from)] <- x[!(x %in% from)]
  res
}

# Recursively flatten a nested list of data frames into a flat list
flattenlist <- function(x) {
  morelists <- sapply(x, function(xprime) class(xprime)[1] == "list")
  out <- c(x[!morelists], unlist(x[morelists], recursive = FALSE))
  if (sum(morelists)) {
    Recall(out)
  } else {
    return(out)
  }
}

# Last observation carried forward
locf <- function(x) {
  good <- !is.na(x)
  positions <- seq(length(x))
  good.positions <- good * positions
  last.good.position <- cummax(good.positions)
  last.good.position[last.good.position == 0] <- NA
  x[last.good.position]
}


# ----------------------------------------------------------------------------
# Builders for the NONMEM event records.
#
# The original code repeated each of these data.frame() constructions many
# times (read-dose vs. argument-driven, dosing history vs. steady state, mg vs.
# mg/kg). They differ only in which subject vector / repeat counts are supplied,
# so they are factored here. "sid"/"ns" are the subject-id vector and the number
# of subjects to recycle over for the regimen being processed.
# ----------------------------------------------------------------------------

# CMT column for a dosing record: primary compartment, interleaved with a second
# compartment when sec.dose.cmt is supplied. "reps" is how many doses each
# compartment value spans (num.dose, num.dose.ss, or 1 for an ADDL record).
build.cmt <- function(dose.cmt.i, sec.dose.cmt.i, reps) {
  c(rbind(rep(dose.cmt.i, each = reps),
          if (all(is.na(sec.dose.cmt.i))) NULL else rep(sec.dose.cmt.i, each = reps)))
}

# RATE column for a dosing record. A negative inf.time is passed straight through
# to RATE; otherwise RATE = dose / inf.time. A second compartment's rate (e.g. a
# duration flag) is interleaved when sec.inf.time is supplied.
build.rate <- function(inf.time.i, sec.inf.time.i, dose.i, reps) {
  c(rbind(rep(if (all(na.omit(inf.time.i < 0))) inf.time.i else dose.i / inf.time.i, each = reps),
          if (all(is.na(sec.inf.time.i))) NULL else rep(sec.inf.time.i, each = reps)))
}

# Observation (concentration) records. "doseno" is the DOSENO vector (already
# expanded to one value per observation row before subject recycling).
build.conc <- function(sid, ns, st.time.i, n.st.i, obs.cmt, n.obs.cmt, di.i, doseno) {
  data.frame(ID     = rep(sid, each = n.st.i * n.obs.cmt),
             TIME   = rep(rep(st.time.i, each = n.obs.cmt), ns),
             DV     = NA,
             EVID   = 0,
             CMT    = rep(rep(obs.cmt, n.st.i), ns),
             DOSENO = rep(doseno, ns),
             EVENT  = 1,
             DI     = di.i)
}

# Dosing-history records. "amt" is supplied by the caller because it differs
# between argument-driven dosing (a scalar dose, recycled) and read-dose-info
# (a per-subject dose vector).
build.dh <- function(sid, ns, dose.time.i, n.dose.cmt, num.dose.i, pre.num.dose, di.i,
                     amt, dose.cmt.i, sec.dose.cmt.i, inf.time.i, sec.inf.time.i, dose.i) {
  data.frame(ID     = rep(sid, each = num.dose.i * n.dose.cmt),
             TIME   = rep(rep(dose.time.i, each = n.dose.cmt), ns),
             EVID   = 1,
             CMT    = build.cmt(dose.cmt.i, sec.dose.cmt.i, num.dose.i),
             DOSENO = rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose.i), each = n.dose.cmt), ns),
             EVENT  = 4,
             DI     = di.i,
             AMT    = amt,
             RATE   = build.rate(inf.time.i, sec.inf.time.i, dose.i, num.dose.i))
}

# AUC control points (compartment off then on at each dose) for the dosing history.
build.cmt.switch <- function(sid, ns, dose.time.i, num.dose.i, pre.num.dose, di.i, auc.cmt) {
  num.auc.cmt <- length(auc.cmt)
  data.frame(ID     = rep(sid, each = num.dose.i * num.auc.cmt * 2),
             TIME   = rep(rep(dose.time.i, each = num.auc.cmt * 2), ns),
             CMT    = rep(rep(c(-auc.cmt, auc.cmt), num.dose.i), ns),
             EVID   = 2,
             DOSENO = rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose.i), each = num.auc.cmt * 2), ns),
             EVENT  = rep(rep(c(rep(2, num.auc.cmt), rep(3, num.auc.cmt)), num.dose.i), ns),
             DI     = rep(rep(di.i, num.dose.i * num.auc.cmt * 2), ns))
}

# AUC control points at steady state (DOSENO is a single value; doses span addl..addl+num.dose.ss).
build.cmt.switch.ss <- function(sid, ns, num.dose.ss.i, addl, di.i, auc.cmt, doseno) {
  num.auc.cmt <- length(auc.cmt)
  data.frame(ID     = rep(sid, each = num.dose.ss.i * num.auc.cmt * 2),
             TIME   = rep(rep(seq(addl * di.i, di.i * (num.dose.ss.i + addl - 1), di.i), each = num.auc.cmt * 2), ns),
             CMT    = rep(rep(c(-auc.cmt, auc.cmt), num.dose.ss.i), ns),
             EVID   = 2,
             DOSENO = doseno,
             EVENT  = rep(rep(c(rep(2, num.auc.cmt), rep(3, num.auc.cmt)), num.dose.ss.i), ns),
             DI     = di.i)
}

# Merge the event table ("coredf" = concentration + dosing-history records) with the
# covariates. Handles both the one-record-per-ID case and the time-varying case (where
# the merge is done in two steps and gaps are back-filled by EVID=2 / locf).
merge.covariates <- function(coredf, cov.df, cov.df.long, tv.cov, remove.col, id.var, time.var) {
  if (is.null(cov.df.long)) {
    if (length(remove.col) > 0) out <- merge(coredf, cov.df[-remove.col], by.x = 'ID', by.y = id.var, all = T)
    else                        out <- merge(coredf, cov.df,              by.x = 'ID', by.y = id.var, all = T)
  } else {
    # merge time-varying covariates first (keyed on both ID and TIME)
    if (length(remove.col) > 0) out <- merge(coredf, cov.df.long[-remove.col][c(id.var, time.var, tv.cov)], by.x = c('ID', 'TIME'), by.y = c(id.var, time.var), all = T)
    else                        out <- merge(coredf, cov.df.long[c(id.var, time.var, tv.cov)],              by.x = c('ID', 'TIME'), by.y = c(id.var, time.var), all = T)

    # then the non-time-varying covariates, if any
    if (length(cov.df) > length(cov.df[c(id.var, time.var, tv.cov)])) {
      tv.col        <- na.omit(match(c(time.var, tv.cov), names(cov.df)))
      remove.tv.col <- na.omit(c(remove.col, tv.col))
      if (length(remove.tv.col) > 0) out <- merge(out, cov.df[-remove.tv.col], by.x = 'ID', by.y = id.var, all = T)
      else                           out <- merge(out, cov.df,                 by.x = 'ID', by.y = id.var, all = T)
    }

    # fill the blanks created by the outer merge
    out[is.na(out$EVID), 'EVID'] <- 2
    col.with.blanks <- c('CMT', 'DOSENO', 'EVENT', 'DI', tv.cov)
    for (cb in col.with.blanks) out[[cb]] <- locf(out[[cb]])
  }
  out
}


# ----------------------------------------------------------------------------
# Main function
# ----------------------------------------------------------------------------
create.sim.dataset <- function(read.dose.info, di, num.dose, dose, dose.unit='mg', dose.cmt=1, sec.dose.cmt=NA, obs.cmt=1, smp.time, smp.time.as.is=F,
                               inf.time=NA, sec.inf.time=NA, auc.interval=TRUE, auc.cmt=3, cov.df, id.var='ID', bw.var='BBWT', time.var = NULL, tv.bw.var=NULL, bw.change.threshold = NULL,
                               output.suffix, export=TRUE, export.path='../../nm/0data/',
                               ss=TRUE, ss.with.addl=FALSE, addl=NULL, auc.ss.interval=FALSE, only.ss=FALSE, smp.time.ss=NULL, num.dose.ss=NULL, skip.dir.check=FALSE, quote=F){

  #check if the export directory exists, if not, create it.
  if(skip.dir.check==FALSE){
    if (!dir.exists(export.path)) {
      cat('\nThe "export.path" does not exist. \nCreating ', export.path,'\n', sep='')
      dir.create(export.path,recursive = TRUE)
    }
  }

  #replace inf.time 0 with NA
  inf.time[inf.time==0] <- NA
  sec.inf.time[sec.inf.time==0] <- NA

  #where a sec.inf.time is not NA, the corresponding position of sec.dose.cmt cannot be NA either
  if(any(is.na(sec.dose.cmt[!is.na(sec.inf.time)]))) stop('When an element of sec.inf.time is not NA, the corresponding position sec.dose.cmt cannot be NA either')

  if(!is.data.frame(cov.df)) stop('"cov.df" must be a data frame')
  #remove data.table class if cov.df is created by data.table package. The current code does not work with data.table
  if(is.data.table(cov.df)) class(cov.df) <- 'data.frame'

  #check if there is any NA or -99 in the cov.df
  if(sum(is.na(cov.df), na.rm=T) > 0) message(paste(paste(names(which(colSums(is.na(cov.df)) > 0)), collapse =' '), 'contains NA', sep =' '))
  if(sum(cov.df==-99, na.rm=T) > 0) message(paste(paste(names(which(colSums(cov.df==-99) > 0)), collapse =' '), 'contains -99', sep =' '))

  if(is.null(tv.bw.var) & !is.null(bw.change.threshold)) stop('"tv.bw.var" cannot be empty when "bw.change.threshold" is not NULL')

  # ---- Handle a cov.df with multiple records per ID (time-varying covariates) ----
  if(length(unique(cov.df[[id.var]])) != dim(cov.df)[1]) {
    message('"cov.df" contains multiple records per ID')

    if(is.null(time.var)) {
      time.var <- 'TIME'
      message('"time.var" is defaulted to "TIME"')
    }

    #check if the variables really exist in the cov.df
    if(is.numeric(bw.change.threshold)) names.check <- c(id.var, bw.var, time.var, tv.bw.var)
    else names.check <- c(id.var, bw.var, time.var, tv.bw.var, bw.change.threshold)
    names.check <- names.check[!sapply(names.check, is.null)] #get the variables when they are provided
    names.no.match <- names.check[is.na(match(names.check, names(cov.df)))] #extract the variables that don't match any names in the cov.df
    if(length(names.no.match)>0) stop(paste(paste(names.no.match, collapse = ' '), 'do not exist in "cov.df"', sep = ' '))

    #if bw.change.threshold argument is set, add a column for purpose of adjusting dose based on the threshold of body weight change
    #compare with previous measurement not the baseline measurement
    #for time varying covariate (e.g. body weight), adjust the dose based on the body weight of the dosing time
    if(!is.null(tv.bw.var)) { #if tv.bw.var is defined, dose adjustment based on tv.bw.var is executed

      #sorting the cov.df in case it is not sorted already
      cov.df <- cov.df[with(cov.df, order(get(id.var), get(time.var))), ]

      #check if the first record of each ID starts with TIME zero.
      unique.cov.df <- cov.df[!duplicated(cov.df[[id.var]]) & cov.df[[time.var]]>0,]
      #If not, force the time of the first record to be zero
      if(dim(unique.cov.df)[1] > 0) {
        cov.df[!duplicated(cov.df[[id.var]]) & cov.df[[time.var]]>0, time.var] <- 0
        message(paste(id.var, paste(unique.cov.df[[id.var]], collapse = ','), 'do not start with time zero. The time of the first record of these subjects is then forced to be zero', sep = ' '))
      }

      if(!is.null(bw.change.threshold)) { #if bw.change.threshold is defined, the dose adjustment will be subject to a threshold change
        #calculate a new variable .DCBW (dose change body weight)
        cov.df$.DCBW <- NA #create a new variable .DCBW to hold the body weight that used for dose adjustment, then update them through the for loop
        anchor.bw <- cov.df$.DCBW[1] <- cov.df[[bw.var]][1] # set the anchor BW to be the first BW value of the first subject
        for (i in 2:dim(cov.df)[1]) {
          if(cov.df[[time.var]][i]!=0){ #
            if(abs(anchor.bw - cov.df[[tv.bw.var]][i])/anchor.bw > if(is.numeric(bw.change.threshold)) bw.change.threshold else cov.df[[bw.change.threshold]][i]) { #BW change greater than the bw.change.threshold
              anchor.bw <- cov.df$.DCBW[i] <- cov.df[[tv.bw.var]][i] #assign the corresponding tv.bw.var to .DCBW and anchor.bw
            }
            else cov.df$.DCBW[i] <- anchor.bw
          }
          else  anchor.bw <- cov.df$.DCBW[i] <- cov.df[[bw.var]][i] # when TIME = 0, set the anchor.bw
        }
      }
    }

    #detect time varying covariates
    #count the number of unique values by each ID for each variable in the cov.df
    tmp1 <- t(sapply(split(cov.df, cov.df[id.var]), function(x){
      sapply(x, function(y) length(na.omit(unique(y)))) #count the number of unique values by id.var
    }))
    #tmp1 is a matrix, check if any column has value greater than 1 in any ID
    tv.cov <- names(which(apply(tmp1, MARGIN=2, function(x) any(x>1))))
    tv.cov <- tv.cov[-match(time.var, tv.cov)] #remove time variable from tv.cov

    cov.df.long <- cov.df #save the original cov.df as the long form
    cov.df <- cov.df[!duplicated(cov.df[[id.var]]), ] # create a data frame with one record per ID

    cat('Time varying covariates:', tv.cov, '\n')
  }
  else {
    cov.df.long <- NULL
    tv.cov <- NULL
  }

  #check if missing values are present in di, dose or dose.unit
  lapply(c(di, dose, dose.unit), function(x) if(!is.na(match(x, names(cov.df))) && any(is.na(cov.df[[x]]))) stop (paste0('"',x,'"', ' cannot contain NA', sep='')))

  #get the size of di vector
  if(any(!is.na(match(di, names(cov.df))))) check.di <- unique(cov.df[[di]])[order(unique(cov.df[[di]]))]
  else check.di <- unlist(di)

  # ---- Validate smp.time / smp.time.ss against the dosing intervals ----
  if(!is.null(smp.time) && !all(is.na(smp.time))) {
    if(!is.list(smp.time)) stop('"smp.time" must be a list where each element is a sampling vector for a unique dosing regimen')

    #first check if the size of "smp.time" match with "di"
    if(length(check.di) != length(smp.time)) stop(paste0('The data frame "cov.df" has ', length(check.di), ' dose interval(s), when "smp.time" has ', length(smp.time), ' level(s).', sep=''))

    if(smp.time.as.is==FALSE) {
	  #reorder to make sure the smp.time vector in ascending order by trough sampling time when read.dose.info=TRUE
	  if(read.dose.info) smp.time <- smp.time[order(sapply(smp.time, max))]

      #check if the "smp.time" includes the trough time
      check.smp.time <- smp.time[!is.na(smp.time)]
      check.di <- check.di[!is.na(smp.time)]
      lapply(1:length(check.di), function(x) if (is.na(match(check.di[[x]], check.smp.time[[x]]))) warning(paste0('sampling vector ', paste0(check.smp.time[[x]],collapse=','), ' does not contain the trough time: ', check.di[[x]], sep='')))

      #check if the max "smp.time" is larger than the respective dosing interval
      lapply(1:length(check.di), function(x) if (max(check.smp.time[[x]], na.rm=T) > check.di[[x]]) warning(paste0('sampling vector ', paste0(check.smp.time[[x]],collapse=','), ' contains the sampling time longer than dosing interval: ', check.di[[x]], sep='')))
    }
  }

  if(!is.null(smp.time.ss) && !all(is.na(smp.time.ss))) {
    if(!is.list(smp.time.ss)) stop('"smp.time.ss" must be a list')

    #first check if the size of "smp.time.ss" match with "di"
    #for smp.time.ss, only need to check when read.dose.info=T
    if(read.dose.info) {
      if(length(check.di) != length(smp.time.ss)) stop(paste0('The data frame "cov.df" has ', length(check.di), ' dose interval(s), when "smp.time.ss" has ', length(smp.time.ss), ' level(s).', sep=''))
    }
    else {
      if(length(smp.time.ss)>1) stop('When "read.dose.info"=FALSE, only one time vector should be provided to "smp.time.ss"')
    }

    if(smp.time.as.is==FALSE) {
	  #reorder to make sure the smp.time vector in ascending order by trough sampling time when read.dose.info=TRUE
      if(read.dose.info) smp.time.ss <- smp.time.ss[order(sapply(smp.time.ss, max))]

      #check if the "smp.time.ss" includes the trough time
      if(read.dose.info){
        check.smp.time.ss <- smp.time.ss[!is.na(smp.time.ss)]
        check.di <- check.di[!is.na(smp.time.ss)]
        lapply(1:length(check.di), function(x) if (is.na(match(check.di[[x]], check.smp.time.ss[[x]]))) warning(paste0('sampling vector ', paste0(check.smp.time.ss[[x]],collapse=','), ' does not contain the trough time: ', check.di[[x]], sep='')))
      }
      else{
        if(is.na(match(di[[length(di)]], unlist(smp.time.ss)))) warning(paste0('sampling vector', paste0(unlist(smp.time.ss),collapse=','), ' does not contain the trough time: ', di[[length(di)]], sep=''))
      }

      #check if the max "smp.time.ss" is larger than the respective dosing interval
      if(read.dose.info) lapply(1:length(check.di), function(x) if (max(check.smp.time.ss[[x]], na.rm=T) > check.di[[x]]) warning(paste0('sampling vector ', paste0(check.smp.time.ss[[x]],collapse=','), ' contains the sampling time longer than dosing interval: ', check.di[[x]], sep='')))
      else {if(max(unlist(smp.time.ss), na.rm=T) > di[[length(di)]]) warning(paste0('sampling vector ', paste0(unlist(smp.time.ss),collapse=','), ' contains the sampling time longer than dosing interval: ', di[[length(di)]], sep=''))}
    }
  }

  #check if id.var contains NA
  if(any(is.na(cov.df[[id.var]]))) warning(paste0(id.var, ' in the data frame "cov.df"', ' contains NA', sep=''))

  #check if the reserved names are used in the cov.df
  if(is.null(cov.df.long))  reserved.names <- c('TIME','EVID','CMT','DOSENO','EVENT','DI','DV','AMT','RATE','DAY','WEEK','NDOSESIM','DOSEUNTSIM')
  else reserved.names <- c('EVID','CMT','DOSENO','EVENT','DI','DV','AMT','RATE','DAY','WEEK','NDOSESIM','DOSEUNTSIM')

  if(any(!is.na(match(reserved.names, names(cov.df))))) stop(paste0('reserved names ', reserved.names[!is.na(match(reserved.names, names(cov.df)))], ' are used in the "cov.df"', sep=''))

  #print a warning message if length of obs.cmt > 1
  if(length(obs.cmt) > 1) message(paste0("\n", "THIS IS NOT AN ERROR.", "\n", '"obs.cmt" has ', length(obs.cmt), ' values: ', paste0(obs.cmt, collapse=', '), "\n", 'The values of "obs.cmt" are NOT in pair with the values of other dosing arguments (e.g. "di", "dose", "dose.unit" etc.).', "\n", 'They simply add multiple observation compartments for all sampling time points after each dose.', sep=''))

  # di, dose, dose.unit, inf.time will be redefined later on, save the name for future use
  di.name <- di
  dose.name <- dose
  num.dose.name <- num.dose
  dose.unit.name <- dose.unit
  dose.cmt.name <- dose.cmt
  sec.dose.cmt.name <- sec.dose.cmt
  inf.time.name <- inf.time
  sec.inf.time.name <- sec.inf.time

  # ---- When the dosing arguments are column names, read those values out of cov.df ----
  if(read.dose.info) {
    if(length(di)!=1 || length(dose)!=1 || length(num.dose)!=1 || length(dose.unit)!=1 || length(inf.time)!=1 || length(sec.inf.time)!=1 || length(dose.cmt)!=1 || length(sec.dose.cmt)!=1) stop ('When read.dose.info=TRUE, "di", "dose", "num.dose","dose.unit", "inf.time", "sec.inf.time", "dose.cmt" and "sec.dose.cmt"  must be either a single value or a name of a column in "cov.df"')
    else{
      read.index <- na.omit(match(c(di.name, dose.unit.name, dose.cmt.name, num.dose.name), names(cov.df)))
      read.var.name <- names(cov.df)[read.index]

      if(length(read.var.name) > 0) cov.df.split <- split(x=cov.df, f=cov.df[read.var.name], drop=TRUE)
      else cov.df.split <- list(cov.df) # doing nothing but create a list with cov.df as an element

      #check if the inf.time contain both positive and negative values, if yes, split by inf.time.name (because the processing will be different)
      if(!is.na(match(inf.time.name, names(cov.df)))) {
		  out <- list()
		  out <- lapply(cov.df.split, function(x){
			if(length(unique(na.omit(unique(x[[inf.time.name]])>0)))>1) tmp <- split(x, f=factor(x[[inf.time.name]]>0, exclude=NULL), drop = FALSE) else tmp <- x
			out <- c(out, list(tmp))
		  })
		  cov.df.split <- flattenlist(out)
      }

      #build di vector
      if(!is.na(match(di.name, names(cov.df)))){
        di <- lapply(cov.df.split, function(x)unique(x[[di.name]])) # extract di as unique numeric value
        #build smp.time, smp.time.ss vector
        all.di <- unlist(di)
        unique.di <- unique(all.di)[order(unique(all.di))]
        di.index <- match(all.di, unique.di)
        if(!is.null(smp.time)) smp.time <- smp.time[di.index]
        if(!is.null(smp.time.ss)) smp.time.ss <- smp.time.ss[di.index]
      }
      else di <- di # di will be called as unique value later on

      #build dose vector (NOT a sorting variable)
      if(!is.na(match(dose.name, names(cov.df)))) dose <- lapply(cov.df.split, function(x)x[[dose.name]])
      else dose <- lapply(sapply(cov.df.split, nrow), function(x)rep(dose, x)) # dose will be called as a vector later on

      #build num.dose vector (NOT a sorting variable)
      if(!is.na(match(num.dose.name, names(cov.df)))) num.dose <- lapply(cov.df.split, function(x)unique(x[[num.dose.name]]))
      else num.dose <- num.dose # num.dose will be called as a unique value later on

      #build dose.unit vector
      if(!is.na(match(dose.unit.name, names(cov.df)))) dose.unit <- lapply(cov.df.split, function(x)unique(x[[dose.unit.name]]))
      else dose.unit <- dose.unit # dose unit will be called as unique value later on

      #build dose.cmt vector
      if(!is.na(match(dose.cmt.name, names(cov.df)))) dose.cmt <- lapply(cov.df.split, function(x)x[[dose.cmt.name]])
      else dose.cmt <- lapply(sapply(cov.df.split, nrow), function(x)rep(dose.cmt, x)) # dose.cmt will be called as a vector later on in a length of number of subjects

      #build sec.dose.cmt vector
      if(!is.na(match(sec.dose.cmt.name, names(cov.df)))) sec.dose.cmt <- lapply(cov.df.split, function(x)x[[sec.dose.cmt.name]])
      else sec.dose.cmt <- lapply(sapply(cov.df.split, nrow), function(x)rep(sec.dose.cmt, x)) # sec.dose.cmt will be called as a vector later on

      #build inf.time vector
      if(!is.na(match(inf.time.name, names(cov.df)))) inf.time <- lapply(cov.df.split, function(x)x[[inf.time.name]])
      else inf.time <- lapply(sapply(cov.df.split, nrow), function(x)rep(inf.time, x)) # inf.time will be called as a vector later on

      #build sec.inf.time vector
      if(!is.na(match(sec.inf.time.name, names(cov.df)))) sec.inf.time <- lapply(cov.df.split, function(x)x[[sec.inf.time.name]])
      else sec.inf.time <- lapply(sapply(cov.df.split, nrow), function(x)rep(sec.inf.time, x)) # sec.inf.time will be called as a vector later on

    }
  }

  # ---- Recycle the per-regimen arguments to a common length (max.len) ----
  #get the maximum length of the arguments that may be in vectors
  max.len <- max(sapply(list(di, num.dose, dose, dose.unit, dose.cmt, sec.dose.cmt, smp.time, inf.time, auc.interval), length))
  cat('The number of data frames processed:', max.len, '\n')

  recycle <- function(x, name) {
    if (max.len > 1 && length(x) == 1) x <- rep(x, max.len)
    if (length(x) != max.len) stop(paste0('The length of "', name, '" is not ', max.len, sep=''))
    x
  }
  di           <- recycle(di,           'di')
  num.dose     <- recycle(num.dose,     'num.dose')
  dose         <- recycle(dose,         'dose')
  if(!is.null(smp.time)) smp.time <- recycle(smp.time, 'smp.time')
  dose.cmt     <- recycle(dose.cmt,     'dose.cmt')
  sec.dose.cmt <- recycle(sec.dose.cmt, 'sec.dose.cmt')
  inf.time     <- recycle(inf.time,     'inf.time')
  sec.inf.time <- recycle(sec.inf.time, 'sec.inf.time')
  dose.unit    <- recycle(dose.unit,    'dose.unit')
  auc.interval <- recycle(auc.interval, 'auc.interval')

  #set the number of observation compartment(s)
  n.obs.cmt <- length(obs.cmt)

  #set total number of subjects and unique ID
  if(read.dose.info) {
    n.subj <- sapply(cov.df.split, function(x)dim(x)[1])
    subjID <- lapply(cov.df.split, function(x)x[[id.var]])
    if(length(n.subj)!=max.len) stop(paste0('The length of "n.subj" is not ', max.len, sep=''))
    if(length(subjID)!=max.len) stop(paste0('The length of "subjID" is not ', max.len, sep=''))
  }
  else {
    n.subj <- dim(cov.df)[1]
    subjID <- cov.df[[id.var]]
  }

  # ==========================================================================
  # Construct the dataset for the dosing history (the non-steady-state dataset)
  # ==========================================================================
  if (!only.ss) {

    #define empty lists of vectors used for the contruction of the dataset
    dose.time <- rep(list(vector()), max.len)  # all dosing time points
    st.time <- rep(list(vector()), max.len)    # all sampling time points
    n.st <- rep(list(vector()), max.len)       # total number of sampling time points (the whole time course)
    n.st.di <- rep(list(vector()), max.len)    # the number of sampling time points (within a dosing interval)

    #initialize
    next.regimen.start <- 0     # next regimen start time
    pre.num.dose <- 0           # the accumulative number of doses from the previous regimens
    conc <- data.frame()        # empty data frame for concentration data
    dh <- data.frame()          # empty data frame for dosing history data
    cmt.switch <- data.frame()  # empty data frame for CMT control switch

    for (i in 1:max.len) {

      #subject vector / count for this regimen (a per-split list when reading from cov.df)
      sid <- if(read.dose.info) subjID[[i]] else subjID
      ns  <- if(read.dose.info) n.subj[[i]] else n.subj

      dose.time[[i]] <- seq(next.regimen.start, next.regimen.start + di[[i]] * (num.dose[[i]] - 1), di[[i]])
      if(any(is.na(smp.time[[i]]))) st.time[[i]] <- NA
      else {
        if(smp.time.as.is) st.time[[i]] <- smp.time[[i]]
        else st.time[[i]] <- rep(smp.time[[i]], num.dose[[i]]) + rep(dose.time[[i]], each=length(smp.time[[i]]))
      }
      n.st[[i]] <- length(st.time[[i]])
      n.st.di[[i]] <- length(smp.time[[i]])

      #calculate the number of dosing compartments for each iteration (1 or 2)
      n.dose.cmt <- if(all(is.na(sec.dose.cmt[[i]]))) 1 else 2

      #prepare the concentration data frame
      #Add "EVENT" as a sorting variable to sort different events at the same time
      #EVENT=1 observation, EVENT=2 turn off compartment, EVENT=3 turn on compartment, EVENT=4 dosing
      if(!any(is.na(st.time[[i]]))) {
        #DOSENO assignment for the observation rows
        if(smp.time.as.is)
          doseno <- rep(ifelse(ceiling(st.time[[i]]/di[[i]])==0, 1, ifelse(ceiling((st.time[[i]]-next.regimen.start)/di[[i]])<=num.dose[[i]], ceiling((st.time[[i]]-next.regimen.start)/di[[i]])+pre.num.dose, pre.num.dose+num.dose[[i]])), each=n.obs.cmt)
        else
          doseno <- rep((pre.num.dose + 1):(pre.num.dose + num.dose[[i]]), each=n.st.di[[i]]*n.obs.cmt)
        conc <- rbind(conc, build.conc(sid, ns, st.time[[i]], n.st[[i]], obs.cmt, n.obs.cmt, di[[i]], doseno))
      }

      #construct dosing history
      #AMT is a per-subject vector when reading dose from cov.df, otherwise a single dose recycled
      if(read.dose.info) amt <- rep(rep(dose[[i]], each=n.dose.cmt), each=num.dose[[i]])
      else amt <- dose[[i]]
      dh.i <- build.dh(sid, ns, dose.time[[i]], n.dose.cmt, num.dose[[i]], pre.num.dose, di[[i]],
                       amt, dose.cmt[[i]], sec.dose.cmt[[i]], inf.time[[i]], sec.inf.time[[i]], dose[[i]])

      if (grepl("\\/", ignore.case=TRUE, dose.unit[[i]])) {
        #dose.unit is per body weight (e.g. mg/kg): scale AMT (and RATE) by body weight
        cov.bw <- if(read.dose.info) cov.df.split[[i]][c(id.var, bw.var)] else cov.df[c(id.var, bw.var)]
        dh.i <- merge(dh.i, cov.bw, by.x='ID', by.y=id.var, all.x=T)
        dh.i$AMT <- dh.i[['AMT']] * dh.i[[bw.var]]
        #get the index where RATE is not NA or negative
        rate.index <- !is.na(dh.i$RATE) & dh.i$RATE>0
        if(read.dose.info) rate.denom <- rep(inf.time[[i]], each=num.dose[[i]]*n.dose.cmt)
        else rate.denom <- rep(inf.time[[i]], n.dose.cmt*num.dose[[i]]*n.subj)
        dh.i$RATE[rate.index] <- dh.i$AMT[rate.index] / rate.denom[rate.index]
        dh.i <- dh.i[-match(c(bw.var), names(dh.i))]
      }
      dh <- rbind(dh, dh.i)

      if(auc.interval[[i]]) {
        cmt.switch <- rbind(cmt.switch, build.cmt.switch(sid, ns, dose.time[[i]], num.dose[[i]], pre.num.dose, di[[i]], auc.cmt))
        #remove the control points for the first dose
        cmt.switch <- cmt.switch[cmt.switch$TIME != dose.time[[1]][1], ]
        dh <- merge(dh, cmt.switch, all=T)
      }

      if(read.dose.info){
        next.regimen.start <- 0
        pre.num.dose <- 0
      }
      else{
        next.regimen.start <- next.regimen.start + di[[i]] * num.dose[[i]]
        pre.num.dose <- pre.num.dose + num.dose[[i]]
      }
    }

    #merge all dataset together
    remove.col <- na.omit(match(c(di.name, dose.name, dose.unit.name), names(cov.df))) #get the "di", "dose" and "dose.unit" columns in the cov.df
    sim <- merge.covariates(merge(conc, dh, all=T), cov.df, cov.df.long, tv.cov, remove.col, id.var, time.var)

    sim$DAY <- ceiling(sim$TIME/24)
    sim[sim$DAY==0,]$DAY <- 1
    sim$WEEK <- ceiling(sim$DAY/7)

    #Add NDOSESIM and DOSEUNTSIM
    if(read.dose.info) {
      if(!is.na(match(dose.name, names(cov.df)))) sim$NDOSESIM <- cov.df[[dose.name]][match(sim$ID, cov.df[[id.var]])]  #  help identify an error that cause to produce a NDOSESIM list instead of vector
      else sim$NDOSESIM <- unique(unlist(dose)) # dose, dose.unit may be expanded to a vector with length > 1
      if(!is.na(match(dose.unit.name, names(cov.df)))) sim$DOSEUNTSIM <- cov.df[[dose.unit.name]][match(sim$ID, cov.df[[id.var]])]
      else sim$DOSEUNTSIM <- unique(unlist(dose.unit))
    }
    else{
      sim$NDOSESIM <- map(sim$DOSENO, from=seq_len(sum(num.dose)), to=unlist(lapply(1:length(num.dose), function(x)rep(dose[[x]], num.dose[[x]]))))
      sim$DOSEUNTSIM <- map(sim$DOSENO, from=seq_len(sum(num.dose)), to=unlist(lapply(1:length(num.dose), function(x)rep(dose.unit[[x]], num.dose[[x]]))))
    }

    #dose adjusted by tv.bw.var directly or .DCBW
    if(!is.null(tv.bw.var)) {
      sim[!is.na(sim$AMT), 'AMT'] <- with(sim[!is.na(sim$AMT), ], AMT * (if(!is.null(bw.change.threshold)) .DCBW else get(tv.bw.var)) / get(bw.var))
      if(any(!is.na(sim$RATE))) sim[!is.na(sim$RATE), 'RATE'] <- with(sim[!is.na(sim$RATE), ], RATE * (if(!is.null(bw.change.threshold)) .DCBW else get(tv.bw.var)) /get(bw.var))
    }

    sim$DV <- as.numeric(sim$DV)
    sim <- char.col.at.last(sim)

    #remove RATE if all values are NA
    if(all(is.na(sim$RATE))) sim <- sim[-match('RATE', names(sim))]

    #remove AMT=0 records when EVID=1
    sim <- sim[!(sim$EVID==1 & sim$AMT==0),]

    #reorder
    sim <- with(sim, sim[order(ID,TIME,EVENT,DOSENO),])

    attr(sim, 'class') <- c(class(sim), 'js.csd')

    if (export) {
      write.csv(sim, file=file.path(export.path, paste0(if(is.null(output.suffix)) 'sim' else output.suffix, '.csv')), row.names=F, na='.', quote=quote)
    }

  }

  # ==========================================================================
  # Construct a dataset for steady state
  # ==========================================================================
  if(ss|ss.with.addl) {

    if(is.null(num.dose.ss) || is.na(num.dose.ss)) num.dose.ss <-1
    if(is.character(num.dose.name) && num.dose.ss==num.dose.name) num.dose.ss <- num.dose

    if(max.len>1 && length(num.dose.ss)==1) num.dose.ss <- rep(num.dose.ss, max.len)
    if(length(num.dose)!=max.len) stop(paste0('The length of "num.dose.ss" is not ', max.len, sep=''))

    if(read.dose.info){

      #Get the sampling time points for steady state
      if(is.null(smp.time.ss) || all(is.na(smp.time.ss))) smp.time.ss <- smp.time
      else {
        if(!is.list(smp.time.ss)) stop('"smp.time.ss" must be a list')
        if(length(smp.time.ss) != max.len) stop(paste0('The length of "smp.time.ss" is not ', max.len, sep=''))
      }

      #define empty lists of vectors used for the contruction of the dataset
      dose.time.ss <- rep(list(vector()), max.len)
      st.time.ss <- rep(list(vector()), max.len)
      n.st.ss <- rep(list(vector()), max.len)

      #initialize
      conc.ss <- data.frame()        # empty data frame for concentration data
      dh.ss <- data.frame()          # empty data frame for dosing history data
      cmt.switch.ss <- data.frame()  # empty data frame for CMT control switch at steady state

      for (i in 1:max.len) {
        sid <- subjID[[i]]
        ns  <- n.subj[[i]]

        if(smp.time.as.is) st.time.ss[[i]] <- smp.time.ss[[i]]
        else {
          # when ss.with.addl=TRUE, the st.time.ss also includes the sampling time course up to the num.dose.ss
          if(ss.with.addl) st.time.ss[[i]] <- rep(smp.time.ss[[i]], num.dose.ss[[i]]) + rep(seq(addl*di[[i]], di[[i]]*(num.dose.ss[[i]]+addl-1), di[[i]]), each=length(smp.time.ss[[i]]))
          else st.time.ss[[i]] <- rep(smp.time.ss[[i]], num.dose.ss[[i]]) + rep(seq(0, di[[i]]*(num.dose.ss[[i]]-1), di[[i]]), each=length(smp.time.ss[[i]]))
        }

        n.st.ss[[i]] <-  length(st.time.ss[[i]])

        #calculate the number of dosing compartments for each iteration (1 or 2)
        n.dose.cmt <- if(all(is.na(sec.dose.cmt[[i]]))) 1 else 2

        #DOSENO constant for steady-state records
        doseno.ss <- if(ss.with.addl) addl+1 else 500

        if(ss.with.addl) dose.time.ss[[i]] <- 0
        else dose.time.ss[[i]] <- seq(0, di[[i]]*(num.dose.ss[[i]] - 1), di[[i]])

        #construct concentration data frame at steady state
        if(!any(is.na(st.time.ss[[i]]))) conc.ss <- rbind(conc.ss, build.conc(sid, ns, st.time.ss[[i]], n.st.ss[[i]], obs.cmt, n.obs.cmt, di[[i]], doseno.ss))

        #dosing history for steady state
        # reps spans num.dose.ss doses, except an ADDL record collapses them to a single row (reps=1)
        ss.reps <- if(ss.with.addl) 1 else num.dose.ss[[i]]
        dh.ss.i <- data.frame(ID=rep(rep(sid, each=n.dose.cmt), each=ss.reps), TIME=rep(rep(dose.time.ss[[i]], each=n.dose.cmt), ns), EVID=1,
                              CMT=build.cmt(dose.cmt[[i]], sec.dose.cmt[[i]], ss.reps),
                              DOSENO=doseno.ss, EVENT=4, DI=di[[i]])
        #SS/ADDL/II columns differ between the ADDL and the SS approaches
        if(ss.with.addl) { dh.ss.i$ADDL <- addl+num.dose.ss[[i]]; dh.ss.i$II <- di[[i]] }
        else             { dh.ss.i$SS   <- if(n.dose.cmt==1) 1 else rep(rep(c(1,2), num.dose.ss[[i]]), ns); dh.ss.i$II <- di[[i]] }
        dh.ss.i$AMT  <- rep(rep(dose[[i]], each=n.dose.cmt), each=ss.reps)
        dh.ss.i$RATE <- build.rate(inf.time[[i]], sec.inf.time[[i]], dose[[i]], ss.reps)

        if (grepl("\\/", ignore.case=TRUE, dose.unit[[i]])) {
          dh.ss.i <- merge(dh.ss.i, cov.df.split[[i]][c(id.var, bw.var)], by.x='ID', by.y=id.var, all.x=T)
          dh.ss.i$AMT <- dh.ss.i[['AMT']] * dh.ss.i[[bw.var]]
          #get the index where RATE is not NA or negative
          rate.index <- !is.na(dh.ss.i$RATE) & dh.ss.i$RATE>0
          if(ss.with.addl) rate.denom <- rep(inf.time[[i]], each=n.dose.cmt)
          else rate.denom <- rep(inf.time[[i]], each=num.dose.ss[[i]]*n.dose.cmt)
          dh.ss.i$RATE[rate.index] <- dh.ss.i$AMT[rate.index] / rate.denom[rate.index]
          dh.ss.i <- dh.ss.i[-match(c(bw.var), names(dh.ss.i))]
        }
        dh.ss <- rbind(dh.ss, dh.ss.i)

        if(auc.ss.interval){
          if(ss.with.addl==FALSE) stop('auc.ss.interval==TRUE can only work when ss.with.addl==TRUE')
          cmt.switch.ss <- rbind(cmt.switch.ss, build.cmt.switch.ss(sid, ns, num.dose.ss[[i]], addl, di[[i]], auc.cmt, doseno.ss))
          dh.ss <- merge(dh.ss, cmt.switch.ss, all=T)
        }

      }
    }

    #dosing information is provided through arguments directly
    else{
      regimen.ss <- length(di) #take the last regimen as steady state

      if(is.null(smp.time.ss) || all(is.na(smp.time.ss))) smp.time.ss <- smp.time[[regimen.ss]]
      else smp.time.ss <- unlist(smp.time.ss)

      if(smp.time.as.is) st.time.ss <- smp.time.ss
      else {
        if(ss.with.addl) st.time.ss <- rep(smp.time.ss, num.dose.ss[[regimen.ss]]) + rep(seq(addl*di[[regimen.ss]], di[[regimen.ss]]*(num.dose.ss[[regimen.ss]]+addl-1), di[[regimen.ss]]), each=length(smp.time.ss))
        else st.time.ss <- rep(smp.time.ss, num.dose.ss[[regimen.ss]]) + rep(seq(0, di[[regimen.ss]]*(num.dose.ss[[regimen.ss]]-1), di[[regimen.ss]]), each=length(smp.time.ss))
      }

      n.st.ss <- length(st.time.ss) # the number of sampling time points at steady state

      #calculate the number of dosing compartments for each iteration (1 or 2)
      n.dose.cmt <- if(all(is.na(sec.dose.cmt[[regimen.ss]]))) 1 else 2

      doseno.ss <- if(ss.with.addl) addl+1 else 500

      if(ss.with.addl) dose.time.ss <- 0
      else dose.time.ss <- seq(0, di[[regimen.ss]]*(num.dose.ss[[regimen.ss]] - 1), di[[regimen.ss]])

      #construct the concentration data frame at steady state
      conc.ss <- build.conc(subjID, n.subj, st.time.ss, n.st.ss, obs.cmt, n.obs.cmt, di[[regimen.ss]], doseno.ss)

      #dosing history for steady state
      if (!grepl("\\/", ignore.case=TRUE, dose.unit[[regimen.ss]])) {
        if(ss.with.addl) dh.ss <- data.frame(ID=rep(subjID, each=n.dose.cmt), TIME=rep(rep(dose.time.ss, each=n.dose.cmt), n.subj), EVID=1, CMT=c(rbind(rep(dose.cmt[[regimen.ss]], n.subj), if(all(is.na(sec.dose.cmt[[regimen.ss]]))) NULL else rep(sec.dose.cmt[[regimen.ss]], n.subj))), DOSENO=doseno.ss, EVENT=4, DI=di[[regimen.ss]], ADDL=addl+num.dose.ss[[regimen.ss]], II=di[[regimen.ss]], AMT=dose[[regimen.ss]], RATE=c(rbind(if(all(na.omit(inf.time[[regimen.ss]]<0))) inf.time[[regimen.ss]] else dose[[regimen.ss]]/inf.time[[regimen.ss]], if(all(is.na(sec.inf.time[[regimen.ss]]))) NULL else sec.inf.time[[regimen.ss]])))
        else dh.ss <- data.frame(ID=rep(rep(subjID, each=n.dose.cmt), each=num.dose.ss[[regimen.ss]]), TIME=rep(rep(dose.time.ss, each=n.dose.cmt), n.subj), EVID=1, CMT=c(rbind(rep(dose.cmt[[regimen.ss]], each=num.dose.ss[[regimen.ss]]*n.subj), if(all(is.na(sec.dose.cmt[[regimen.ss]]))) NULL else rep(sec.dose.cmt[[regimen.ss]], each=num.dose.ss[[regimen.ss]]))), DOSENO=doseno.ss, EVENT=4, DI=di[[regimen.ss]], SS=if(n.dose.cmt==1) 1 else rep(rep(c(1,2), num.dose.ss[[regimen.ss]]), n.subj), II=di[[regimen.ss]], AMT=dose[[regimen.ss]], RATE=c(rbind(if(all(na.omit(inf.time[[regimen.ss]]<0))) inf.time[[regimen.ss]] else dose[[regimen.ss]]/inf.time[[regimen.ss]], if(all(is.na(sec.inf.time[[regimen.ss]]))) NULL else sec.inf.time[[regimen.ss]])))
      }
      else {
        if(ss.with.addl) dh.ss <- data.frame(ID=rep(subjID, each=n.dose.cmt), TIME=rep(rep(dose.time.ss, each=n.dose.cmt), n.subj), EVID=1, CMT=c(rbind(rep(dose.cmt[[regimen.ss]], n.subj), if(all(is.na(sec.dose.cmt[[regimen.ss]]))) NULL else rep(sec.dose.cmt[[regimen.ss]], n.subj))), DOSENO=doseno.ss, EVENT=4, DI=di[[regimen.ss]], ADDL=addl+num.dose.ss[[regimen.ss]], II=di[[regimen.ss]], AMT=dose[[regimen.ss]], RATE=c(rbind(if(all(na.omit(inf.time[[regimen.ss]]<0))) inf.time[[regimen.ss]] else dose[[regimen.ss]]/inf.time[[regimen.ss]], if(all(is.na(sec.inf.time[[regimen.ss]]))) NULL else sec.inf.time[[regimen.ss]])))
        else dh.ss <- data.frame(ID=rep(rep(subjID, each=n.dose.cmt), each=num.dose.ss[[regimen.ss]]), TIME=rep(rep(dose.time.ss, each=n.dose.cmt), n.subj), EVID=1, CMT=c(rbind(rep(dose.cmt[[regimen.ss]], each=num.dose.ss[[regimen.ss]]*n.subj), if(all(is.na(sec.dose.cmt[[regimen.ss]]))) NULL else rep(sec.dose.cmt[[regimen.ss]], each=num.dose.ss[[regimen.ss]]))), DOSENO=doseno.ss, EVENT=4, DI=di[[regimen.ss]], SS=if(n.dose.cmt==1) 1 else rep(rep(c(1,2), num.dose.ss[[regimen.ss]]), n.subj), II=di[[regimen.ss]], AMT=dose[[regimen.ss]], RATE=c(rbind(if(all(na.omit(inf.time[[regimen.ss]]<0))) inf.time[[regimen.ss]] else dose[[regimen.ss]]/inf.time[[regimen.ss]], if(all(is.na(sec.inf.time[[regimen.ss]]))) NULL else sec.inf.time[[regimen.ss]])))
        dh.ss <- merge(dh.ss, cov.df[c(id.var, bw.var)], by.x='ID', by.y=id.var, all.x=T)
        dh.ss$AMT <- dh.ss[['AMT']] * dh.ss[[bw.var]]
        #get the index where RATE is not NA or negative
        rate.index <- !is.na(dh.ss$RATE) & dh.ss$RATE>0
        if(ss.with.addl) dh.ss$RATE[rate.index] <- dh.ss$AMT[rate.index] / rep(inf.time[[regimen.ss]], n.dose.cmt*n.subj)[rate.index] #inf.time[[regimen.ss]] is a single value, num.dose.ss doesn't apply for ss.with.addl
        else dh.ss$RATE[rate.index] <- dh.ss$AMT[rate.index] / rep(inf.time[[regimen.ss]], n.dose.cmt*num.dose.ss[[regimen.ss]]*n.subj)[rate.index]
        dh.ss <- dh.ss[-match(c(bw.var), names(dh.ss))]
      }
      if(auc.ss.interval){
        if(ss.with.addl==FALSE) stop('auc.ss.interval==TRUE can only work when ss.with.addl==TRUE')
        cmt.switch.ss <- build.cmt.switch.ss(subjID, n.subj, num.dose.ss[[regimen.ss]], addl, di[[regimen.ss]], auc.cmt, doseno.ss)
        dh.ss <- merge(dh.ss, cmt.switch.ss, all=T)
      }
    }

    #merge all the data frames
    remove.col <- na.omit(match(c(di.name, dose.name, dose.unit.name), names(cov.df))) #get the "di", "dose" and "dose.unit" columns in the cov.df
    sim.ss <- merge.covariates(merge(conc.ss, dh.ss, all=T), cov.df, cov.df.long, tv.cov, remove.col, id.var, time.var)

    sim.ss$DAY <- ceiling(sim.ss$TIME/24)
    sim.ss$WEEK <- ceiling(sim.ss$DAY/7)

    #Add NDOSESIM and DOSEUNTSIM
    if(read.dose.info) {
      if(!is.na(match(dose.name, names(cov.df)))) sim.ss$NDOSESIM <- cov.df[[ dose.name]][match(sim.ss$ID, cov.df[[id.var]])]
      else sim.ss$NDOSESIM <- unique(unlist(dose))
      if(!is.na(match(dose.unit.name, names(cov.df)))) sim.ss$DOSEUNTSIM <- cov.df[[dose.unit.name]][match(sim.ss$ID, cov.df[[id.var]])]
      else sim.ss$DOSEUNTSIM <- unique(unlist(dose.unit))
    }
    else{
      sim.ss$NDOSESIM <- map(sim.ss$DI, from=di[[regimen.ss]], to=dose[[regimen.ss]])
      sim.ss$DOSEUNTSIM <- map(sim.ss$NDOSESIM, from=dose[[regimen.ss]], to=dose.unit[[regimen.ss]])

      sim.ss$NDOSESIM <- dose[[regimen.ss]]
      sim.ss$DOSEUNTSIM <- dose.unit[[regimen.ss]]
    }

    #dose adjusted by tv.bw.var directly or .DCBW
    if(!is.null(tv.bw.var)) {
      sim.ss[!is.na(sim.ss$AMT), 'AMT'] <- with(sim.ss[!is.na(sim.ss$AMT), ], AMT * (if(!is.null(bw.change.threshold)) .DCBW else get(tv.bw.var)) / get(bw.var))
      if(any(!is.na(sim.ss$RATE))) sim.ss[!is.na(sim.ss$RATE), 'RATE'] <- with(sim.ss[!is.na(sim.ss$RATE), ], RATE * (if(!is.null(bw.change.threshold)) .DCBW else get(tv.bw.var)) /get(bw.var))
    }

    sim.ss$DV <- as.numeric(sim.ss$DV)
    sim.ss <- char.col.at.last(sim.ss)

    #reorder
    sim.ss <- with(sim.ss, sim.ss[order(ID,TIME,EVENT,DOSENO),])

    #remove RATE if all values are NA
    if(all(is.na(sim.ss$RATE))) sim.ss <- sim.ss[-match('RATE', names(sim.ss))]

    attr(sim.ss, 'class') <- c(class(sim.ss), 'js.csd')

    if (export) {
      write.csv(sim.ss, file=file.path(export.path, paste0(if(is.null(output.suffix)) 'sim_ss' else paste0('ss_', output.suffix, '.csv'))), row.names=F, na='.', quote=quote)
    }
  }


  if((ss|ss.with.addl) & only.ss) return(invisible(sim.ss))
  if((!ss & !ss.with.addl) & !only.ss) return(invisible(sim))
  if((ss|ss.with.addl) & !only.ss) return(invisible(list(sim=sim, sim.ss=sim.ss)))
}