#**************************************************************************
# 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
#
#			
#
# NOTES: Additional info on code in script that could be useful to reviewer
#
# PATH TO EXAMPLE: 
#
# 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%)
# QC HISTORY:
# DATE       PROGRAMMER            MODIFICATION
# ---------  --------------------  ------------------------------------------
# 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)


#define some functions to be called within the function
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)) {
    match(col.first,names(data)) -> temp
    setdiff(1:dim(data)[2],temp) -> temp2
    data[c(temp,temp2)] -> data
  }
  
  grep(FALSE,sapply(data,is.numeric)) -> char.col
  grep(TRUE,sapply(data,is.numeric)) -> num.col
  
  data[c(num.col,char.col)] -> data
  return(data)
}

#another function
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
}

#another function
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 carry 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]
}

#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')
  
  #check if the cov.df has more than one record per ID
  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)
  
  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 di, dose, dose.unit and/or inf.time values may be provided through 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  
      
    }
  }  

  #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')
  
  if(max.len>1 && length(di)==1) di <- rep(di, max.len)
  if(length(di) != max.len) stop(paste0('The length of "di" is not ', max.len, sep=''))
  
  if(max.len>1 && length(num.dose)==1) num.dose <- rep(num.dose, max.len)
  if(length(num.dose)!=max.len) stop(paste0('The length of "num.dose" is not ', max.len, sep=''))
  
  if(max.len>1 && length(dose)==1) dose <- rep(dose, max.len)
  if(length(dose)!=max.len) stop(paste0('The length of "dose" is not ', max.len, sep=''))
  
  if(!is.null(smp.time)){
    if(max.len>1 && length(smp.time)==1) smp.time <- rep(smp.time, max.len)
    if(length(smp.time)!=max.len) stop(paste0('The length of "smp.time" is not ', max.len, sep=''))
  }
  
  if(max.len>1 && length(dose.cmt)==1) dose.cmt <- rep(dose.cmt, max.len)
  if(length(dose.cmt)!=max.len) stop(paste0('The length of "dose.cmt" is not ', max.len, sep=''))
  
  #if(all(!is.na(sec.dose.cmt))){
    if(max.len>1 && length(sec.dose.cmt)==1) sec.dose.cmt <- rep(sec.dose.cmt, max.len)
    if(length(sec.dose.cmt)!=max.len) stop(paste0('The length of "sec.dose.cmt" is not ', max.len, sep=''))
  #}
  
  if(max.len>1 && length(inf.time)==1) inf.time <- rep(inf.time, max.len)
  if(length(inf.time)!=max.len) stop(paste0('The length of "inf.time" is not ', max.len, sep=''))
  
  #if(all(!is.na(sec.inf.time))) {
    if(max.len>1 && length(sec.inf.time)==1) sec.inf.time <- rep(sec.inf.time, max.len)
    if(length(sec.inf.time)!=max.len) stop(paste0('The length of "sec.inf.time" is not ', max.len, sep=''))
  #}
  
  if(max.len>1 && length(dose.unit)==1) dose.unit <- rep(dose.unit, max.len)
  if(length(dose.unit)!=max.len) stop(paste0('The length of "dose.unit" is not ', max.len, sep=''))
  
  if(max.len>1 && length(auc.interval)==1) auc.interval <- rep(auc.interval, max.len)
  if(length(auc.interval)!=max.len) stop(paste0('The length of "auc.interval" is not ', max.len, sep=''))
  
  
  #set the number of observation compartment(s) 
  n.obs.cmt <- length(obs.cmt)
  
  #create a variable for the number of dosing compartments
  #if(is.na(sec.dose.cmt)) n.dose.cmt <- 1 else n.dose.cmt <- 2
  
  #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 a dosing history
  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
    dh.tmp <- data.frame()      # empty data frame for temp dosing history data
    cmt.switch <- data.frame()  # empty data frame for CMT control switch
    
    for (i in 1:max.len) {
      
      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
      #print(paste('max.len =', max.len))
      #print(sec.dose.cmt[[i]])
      #print(paste('n.dose.cmt =', n.dose.cmt))
      
          #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]]))) {
              if(read.dose.info){
                if(smp.time.as.is) conc <- rbind(conc, data.frame(ID=rep(subjID[[i]], each=n.st[[i]]*n.obs.cmt), TIME=rep(rep(st.time[[i]], each=n.obs.cmt), n.subj[[i]]), DV=NA, EVID=0, CMT=rep(rep(obs.cmt, n.st[[i]]), n.subj[[i]]), DOSENO=rep(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), n.subj[[i]]), EVENT=1, DI=di[[i]]))
                else conc <- rbind(conc, data.frame(ID=rep(subjID[[i]], each=n.st[[i]]*n.obs.cmt), TIME=rep(rep(st.time[[i]], each=n.obs.cmt), n.subj[[i]]), DV=NA, EVID=0, CMT=rep(rep(obs.cmt, n.st[[i]]), n.subj[[i]]), DOSENO=rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose[[i]]),each=n.st.di[[i]]*n.obs.cmt), n.subj[[i]]), EVENT=1, DI=di[[i]]))
              }
              else{
                
                if(smp.time.as.is) conc <- rbind(conc, data.frame(ID=rep(subjID, each=n.st[[i]]*n.obs.cmt), TIME=rep(rep(st.time[[i]], each=n.obs.cmt), n.subj), DV=NA, EVID=0, CMT=rep(rep(obs.cmt, n.st[[i]]), n.subj), DOSENO=rep(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), n.subj), EVENT=1, DI=di[[i]])) 
                else conc <- rbind(conc, data.frame(ID=rep(subjID, each=n.st[[i]]*n.obs.cmt), TIME=rep(rep(st.time[[i]], each=n.obs.cmt), n.subj), DV=NA, EVID=0, CMT=rep(rep(obs.cmt, n.st[[i]]), n.subj), DOSENO=rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose[[i]]),each=n.st.di[[i]]*n.obs.cmt), n.subj), EVENT=1, DI=di[[i]]))  
              }
          }
          #construct dosing history
          if (!grepl("\\/", ignore.case=TRUE, dose.unit[[i]])){
            if(read.dose.info) dh <- rbind(dh, data.frame(ID=rep(rep(subjID[[i]], each=n.dose.cmt), each=num.dose[[i]]), TIME=rep(rep(dose.time[[i]], each=n.dose.cmt), n.subj[[i]]), EVID=1, CMT=c(rbind(rep(dose.cmt[[i]], each=num.dose[[i]]), if(all(is.na(sec.dose.cmt[[i]]))) NULL else rep(sec.dose.cmt[[i]], each=num.dose[[i]]))), DOSENO=rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose[[i]]), each=n.dose.cmt), n.subj[[i]]), EVENT=4, DI=di[[i]], AMT=rep(rep(dose[[i]], each=n.dose.cmt), each=num.dose[[i]]), RATE=c(rbind(rep(if(all(na.omit(inf.time[[i]]<0))) inf.time[[i]] else dose[[i]]/inf.time[[i]], each=num.dose[[i]]), if(all(is.na(sec.inf.time[[i]]))) NULL else rep(sec.inf.time[[i]], each=num.dose[[i]])))))
            else dh <- rbind(dh, data.frame(ID=rep(subjID, each=num.dose[[i]]*n.dose.cmt), TIME=rep(rep(dose.time[[i]], each=n.dose.cmt), n.subj), EVID=1, CMT=c(rbind(rep(dose.cmt[[i]], each=num.dose[[i]]), if(all(is.na(sec.dose.cmt[[i]]))) NULL else rep(sec.dose.cmt[[i]], each=num.dose[[i]]))), DOSENO=rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose[[i]]), each=n.dose.cmt), n.subj), EVENT=4, DI=di[[i]], AMT=dose[[i]], RATE=c(rbind(rep(if(all(na.omit(inf.time[[i]]<0))) inf.time[[i]] else dose[[i]]/inf.time[[i]], each=num.dose[[i]]), if(all(is.na(sec.inf.time[[i]]))) NULL else rep(sec.inf.time[[i]], each=num.dose[[i]])))))
            
          }
          else{
            if(read.dose.info) {
              dh.tmp <- data.frame(ID=rep(rep(subjID[[i]], each=n.dose.cmt), each=num.dose[[i]]), TIME=rep(rep(dose.time[[i]], each=n.dose.cmt), n.subj[[i]]), EVID=1, CMT=c(rbind(rep(dose.cmt[[i]], each=num.dose[[i]]), if(all(is.na(sec.dose.cmt[[i]]))) NULL else rep(sec.dose.cmt[[i]], each=num.dose[[i]]))), DOSENO=rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose[[i]]), each=n.dose.cmt), n.subj[[i]]), EVENT=4, DI=di[[i]], AMT=rep(rep(dose[[i]], each=n.dose.cmt), each=num.dose[[i]]), RATE=c(rbind(rep(if(all(na.omit(inf.time[[i]]<0))) inf.time[[i]] else dose[[i]]/inf.time[[i]], each=num.dose[[i]]), if(all(is.na(sec.inf.time[[i]]))) NULL else rep(sec.inf.time[[i]], each=num.dose[[i]]))))
              dh.tmp <- merge(dh.tmp, cov.df.split[[i]][c(id.var, bw.var)], by.x='ID', by.y=id.var, all.x=T)
              dh.tmp$AMT <- dh.tmp[['AMT']] * dh.tmp[[bw.var]]
              #get the index where RATE is not NA or negative
              rate.index <- !is.na(dh.tmp$RATE) & dh.tmp$RATE>0
              dh.tmp$RATE[rate.index] <- dh.tmp$AMT[rate.index] / rep(inf.time[[i]], each=num.dose[[i]]*n.dose.cmt)[rate.index]
            }
            else {
              dh.tmp <- data.frame(ID=rep(subjID, each=num.dose[[i]]*n.dose.cmt), TIME=rep(rep(dose.time[[i]], each=n.dose.cmt), n.subj), EVID=1, CMT=c(rbind(rep(dose.cmt[[i]], each=num.dose[[i]]), if(all(is.na(sec.dose.cmt[[i]]))) NULL else rep(sec.dose.cmt[[i]], each=num.dose[[i]]))), DOSENO=rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose[[i]]), each=n.dose.cmt), n.subj), EVENT=4, DI=di[[i]], AMT=dose[[i]], RATE=c(rbind(rep(if(all(na.omit(inf.time[[i]]<0))) inf.time[[i]] else dose[[i]]/inf.time[[i]], each=num.dose[[i]]), if(all(is.na(sec.inf.time[[i]]))) NULL else rep(sec.inf.time[[i]], each=num.dose[[i]]))))
              dh.tmp <- merge(dh.tmp, cov.df[c(id.var, bw.var)], by.x='ID', by.y=id.var, all.x=T)
              dh.tmp$AMT <- dh.tmp[['AMT']] * dh.tmp[[bw.var]]
              #get the index where RATE is not NA or negative
              rate.index <- !is.na(dh.tmp$RATE) & dh.tmp$RATE>0
              dh.tmp$RATE[rate.index] <- dh.tmp$AMT[rate.index] / rep(inf.time[[i]], n.dose.cmt*num.dose[[i]]*n.subj)[rate.index]
            }
            dh.tmp[-match(c(bw.var), names(dh.tmp))] -> dh.tmp
            dh <- rbind(dh, dh.tmp)
          }

          if(auc.interval[[i]]) {
            num.auc.cmt <- length(auc.cmt)
            if(read.dose.info) cmt.switch <- rbind(cmt.switch, data.frame(ID=rep(subjID[[i]], each=num.dose[[i]]*num.auc.cmt*2), TIME=rep(rep(dose.time[[i]], each=num.auc.cmt*2), n.subj[[i]]), CMT=rep(rep(c(-auc.cmt, auc.cmt), num.dose[[i]]), n.subj[[i]]), EVID=2, DOSENO=rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose[[i]]),each=num.auc.cmt*2), n.subj[[i]]), EVENT=rep(rep(c(rep(2, num.auc.cmt),rep(3, num.auc.cmt)), num.dose[[i]]), n.subj[[i]]), DI=rep(rep(di[[i]], num.dose[[i]]*num.auc.cmt*2), n.subj[[i]])))
            else cmt.switch <- rbind(cmt.switch, data.frame(ID=rep(subjID, each=num.dose[[i]]*num.auc.cmt*2), TIME=rep(rep(dose.time[[i]], each=num.auc.cmt*2), n.subj), CMT=rep(rep(c(-auc.cmt, auc.cmt), num.dose[[i]]), n.subj), EVID=2, DOSENO=rep(rep((pre.num.dose + 1):(pre.num.dose + num.dose[[i]]),each=num.auc.cmt*2), n.subj), EVENT=rep(rep(c(rep(2, num.auc.cmt),rep(3, num.auc.cmt)), num.dose[[i]]), n.subj), DI=rep(rep(di[[i]], num.dose[[i]]*num.auc.cmt*2), n.subj)))
            #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
    if(is.null(cov.df.long)) {
        if (length(remove.col) > 0) sim <- merge(merge(conc, dh, all=T), cov.df[-remove.col], by.x='ID', by.y=id.var, all=T)
        else sim <- merge(merge(conc, dh, all=T), cov.df, by.x='ID', by.y=id.var, all=T)
    }
    else{ #when time varying covariates involved, the merge should be done in two steps: tv and non-tv
        #merge tv covariates first, merge only the tv.cov
        if (length(remove.col) > 0) sim <- merge(merge(conc, dh, all=T), 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 sim <- merge(merge(conc, dh, all=T), cov.df.long[c(id.var,time.var,tv.cov)], by.x=c('ID','TIME'), by.y=c(id.var,time.var), all=T)
      
               
        #check if there is any non time varying cov,if yes do the second merge
        if(length(cov.df) > length(cov.df[c(id.var,time.var,tv.cov)])) {
          #get the position of time varying covariates (and time variable) for exclusion
          tv.col <- na.omit(match(c(time.var, tv.cov), names(cov.df))) #tv.cov has been defined at the beginning of the function
          remove.tv.col <- na.omit(c(remove.col, tv.col))
          if (length(remove.tv.col) > 0) sim <- merge(sim, cov.df[-remove.tv.col], by.x='ID', by.y=id.var, all=T)
          else sim <- merge(sim, cov.df, by.x='ID', by.y=id.var, all=T)
        }
        
        
        #fill EVID=2 for the blanks
        sim[is.na(sim$EVID), 'EVID'] <- 2
        #locf to fill in the blanks
        col.with.blanks <- c('CMT','DOSENO','EVENT','DI', tv.cov)
        for(i in 1:length(col.with.blanks)) sim[[col.with.blanks[i]]] <- locf(sim[[col.with.blanks[i]]]) 
    }
    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)
    char.col.at.last(sim) -> 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) || 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
      dh.tmp.ss <- data.frame()      # empty data frame for temp dosing history data
      cmt.switch.ss <- data.frame()  # empty data frame for CMT control switch at steady state
      
      for (i in 1:max.len) {
        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
        
        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, data.frame(ID=rep(subjID[[i]], each=n.st.ss[[i]]*n.obs.cmt), TIME=rep(rep(st.time.ss[[i]], each=n.obs.cmt), n.subj[[i]]), DV=NA, EVID=0, CMT=rep(rep(obs.cmt, n.st.ss[[i]]), n.subj[[i]]), DOSENO=if(ss.with.addl) addl+1 else 500, EVENT=1, DI=di[[i]])) 
        
        #dosing history for steady state
        if (!grepl("\\/", ignore.case=TRUE, dose.unit[[i]])) {
          if(ss.with.addl) dh.ss <- rbind(dh.ss, data.frame(ID=rep(subjID[[i]], each=n.dose.cmt), TIME=rep(rep(dose.time.ss[[i]], each=n.dose.cmt), n.subj[[i]]), EVID=1, CMT=c(rbind(dose.cmt[[i]], if(all(is.na(sec.dose.cmt[[i]]))) NULL else sec.dose.cmt[[i]])), DOSENO=if(ss.with.addl) addl+1 else 500, EVENT=4, DI=di[[i]], ADDL=addl+num.dose.ss[[i]], II=di[[i]], AMT=rep(dose[[i]], each=n.dose.cmt), RATE=c(rbind(if(all(na.omit(inf.time[[i]]<0))) inf.time[[i]] else dose[[i]]/inf.time[[i]], if(all(is.na(sec.inf.time[[i]]))) NULL else sec.inf.time[[i]]))))
          else  dh.ss <- rbind(dh.ss, data.frame(ID=rep(rep(subjID[[i]], each=n.dose.cmt), each=num.dose.ss[[i]]), TIME=rep(rep(dose.time.ss[[i]], each=n.dose.cmt), n.subj[[i]]), EVID=1, CMT=c(rbind(rep(dose.cmt[[i]], each=num.dose.ss[[i]]), if(all(is.na(sec.dose.cmt[[i]]))) NULL else rep(sec.dose.cmt[[i]], each=num.dose.ss[[i]]))), DOSENO=if(ss.with.addl) addl+1 else 500, EVENT=4, DI=di[[i]], SS=if(n.dose.cmt==1) 1 else rep(rep(c(1,2), num.dose.ss[[i]]), n.subj[[i]]), II=di[[i]], AMT=rep(rep(dose[[i]], each=n.dose.cmt), each=num.dose.ss[[i]]), RATE=c(rbind(rep(if(all(na.omit(inf.time[[i]]<0))) inf.time[[i]] else dose[[i]]/inf.time[[i]], each=num.dose.ss[[i]]), if(all(is.na(sec.inf.time[[i]]))) NULL else rep(sec.inf.time[[i]], each=num.dose.ss[[i]])))))
        }
        else{
          if(ss.with.addl) dh.tmp.ss <- data.frame(ID=rep(subjID[[i]], each=n.dose.cmt), TIME=rep(rep(dose.time.ss[[i]], each=n.dose.cmt), n.subj[[i]]), EVID=1, CMT=c(rbind(dose.cmt[[i]], if(all(is.na(sec.dose.cmt[[i]]))) NULL else sec.dose.cmt[[i]])), DOSENO=if(ss.with.addl) addl+1 else 500, EVENT=4, DI=di[[i]], ADDL=addl+num.dose.ss[[i]], II=di[[i]], AMT=rep(dose[[i]], each=n.dose.cmt), RATE=c(rbind(if(all(na.omit(inf.time[[i]]<0))) inf.time[[i]] else dose[[i]]/inf.time[[i]], if(all(is.na(sec.inf.time[[i]]))) NULL else sec.inf.time[[i]])))
          else dh.tmp.ss <- data.frame(ID=rep(rep(subjID[[i]], each=n.dose.cmt), each=num.dose.ss[[i]]), TIME=rep(rep(dose.time.ss[[i]], each=n.dose.cmt), n.subj[[i]]), EVID=1, CMT=c(rbind(rep(dose.cmt[[i]], each=num.dose.ss[[i]]), if(all(is.na(sec.dose.cmt[[i]]))) NULL else rep(sec.dose.cmt[[i]], each=num.dose.ss[[i]]))), DOSENO=if(ss.with.addl) addl+1 else 500, EVENT=4, DI=di[[i]], SS=if(n.dose.cmt==1) 1 else rep(rep(c(1,2), num.dose.ss[[i]]), n.subj[[i]]), II=di[[i]], AMT=rep(rep(dose[[i]], each=n.dose.cmt), each=num.dose.ss[[i]]), RATE=c(rbind(rep(if(all(na.omit(inf.time[[i]]<0))) inf.time[[i]] else dose[[i]]/inf.time[[i]], each=num.dose.ss[[i]]), if(all(is.na(sec.inf.time[[i]]))) NULL else rep(sec.inf.time[[i]], each=num.dose.ss[[i]]))))
          dh.tmp.ss <- merge(dh.tmp.ss, cov.df.split[[i]][c(id.var, bw.var)], by.x='ID', by.y=id.var, all.x=T)
          dh.tmp.ss$AMT <- dh.tmp.ss[['AMT']] * dh.tmp.ss[[bw.var]]
          #get the index where RATE is not NA or negative
          rate.index <- !is.na(dh.tmp.ss$RATE) & dh.tmp.ss$RATE>0
          if(ss.with.addl) dh.tmp.ss$RATE[rate.index] <- dh.tmp.ss$AMT[rate.index] / rep(inf.time[[i]], each=n.dose.cmt)[rate.index]
          else dh.tmp.ss$RATE[rate.index] <- dh.tmp.ss$AMT[rate.index] / rep(inf.time[[i]], each=num.dose.ss[[i]]*n.dose.cmt)[rate.index]
          dh.tmp.ss[-match(c(bw.var), names(dh.tmp.ss))] -> dh.tmp.ss
          dh.ss <- rbind(dh.ss, dh.tmp.ss)
        }
        
        if(auc.ss.interval){
          num.auc.cmt <- length(auc.cmt)
          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, data.frame(ID=rep(subjID[[i]], 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), n.subj[[i]]), CMT=rep(rep(c(-auc.cmt, auc.cmt), num.dose.ss[[i]]), n.subj[[i]]), EVID=2, DOSENO=if(ss.with.addl) addl+1 else 500, EVENT=rep(rep(c(rep(2, num.auc.cmt),rep(3, num.auc.cmt)), num.dose.ss[[i]]), n.subj[[i]]), DI=di[[i]]))
          #cmt.switch.ss <- cmt.switch.ss[cmt.switch.ss$TIME != dose.time.ss[[1]][1], ]
          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
      
      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 <- data.frame(ID=rep(subjID, each=n.st.ss*n.obs.cmt), TIME=rep(rep(st.time.ss, each=n.obs.cmt), n.subj), DV=NA, EVID=0, CMT=rep(rep(obs.cmt, n.st.ss), n.subj), DOSENO=if(ss.with.addl) addl+1 else 500, EVENT=1, DI=di[[regimen.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=if(ss.with.addl) addl+1 else 500, 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=if(ss.with.addl) addl+1 else 500, 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=if(ss.with.addl) addl+1 else 500, 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=if(ss.with.addl) addl+1 else 500, 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[-match(c(bw.var), names(dh.ss))] -> dh.ss
      }
      if(auc.ss.interval){
        num.auc.cmt <- length(auc.cmt)
        if(ss.with.addl==FALSE) stop('auc.ss.interval==TRUE can only work when ss.with.addl==TRUE')
        cmt.switch.ss <- data.frame(ID=rep(subjID, each=num.dose.ss[[regimen.ss]]*num.auc.cmt*2), TIME=rep(rep(seq(addl*di[[regimen.ss]], di[[regimen.ss]]*(num.dose.ss[[regimen.ss]]+addl-1), di[[regimen.ss]]), each=num.auc.cmt*2), n.subj), CMT=rep(rep(c(-auc.cmt, auc.cmt), num.dose.ss[[regimen.ss]]), n.subj), EVID=2, DOSENO=if(ss.with.addl) addl+1 else 500, EVENT=rep(rep(c(rep(2, num.auc.cmt),rep(3, num.auc.cmt)), num.dose.ss[[regimen.ss]]), n.subj), DI=di[[regimen.ss]])
        #cmt.switch.ss <- cmt.switch.ss[cmt.switch.ss$TIME != dose.time.ss[[1]][1], ]
        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
    if(is.null(cov.df.long)){
        if (length(remove.col) > 0) sim.ss <- merge(merge(conc.ss, dh.ss, all=T), cov.df[-remove.col], by.x='ID', by.y=id.var, all=T)
        else sim.ss <- merge(merge(conc.ss, dh.ss, all=T), cov.df, by.x='ID', by.y=id.var, all=T)
    }
    else {#when time varying covariates involved, the merge should be done in two steps: tv and non-tv
        #merge tv covariates first, merge only the tv.cov
        if (length(remove.col) > 0) sim.ss <- merge(merge(conc.ss, dh.ss, all=T), 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 sim.ss <- merge(merge(conc.ss, dh.ss, all=T), cov.df.long[c(id.var,time.var,tv.cov)], by.x=c('ID','TIME'), by.y=c(id.var, time.var), all=T)
    
        #check if there is any non time varying cov,if yes do the second merge
        if(length(cov.df) > length(cov.df[c(id.var,time.var,tv.cov)])) {
          #get the position of time varying covariates (and time variable) for exclusion
          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) sim.ss <- merge(sim.ss, cov.df[-remove.tv.col], by.x='ID', by.y=id.var, all=T)
          else sim.ss <- merge(sim.ss, cov.df, by.x='ID', by.y=id.var, all=T)
        }
        
        #fill EVID=2 for the blanks
        sim.ss[is.na(sim.ss$EVID), 'EVID'] <- 2
        #locf to fill in the blanks
        col.with.blanks <- c('CMT','DOSENO','EVENT','DI', tv.cov)
        for(i in 1:length(col.with.blanks)) sim.ss[[col.with.blanks[i]]] <- locf(sim.ss[[col.with.blanks[i]]]) 
        
    }
    sim.ss$DAY <- ceiling(sim.ss$TIME/24)
    #sim.ss[sim.ss$DAY==0,]$DAY <- 1
    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)
    char.col.at.last(sim.ss) -> 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)))
}