#Simulating from a fixed network

[yiberkeley]

The current simulation seems to only allow us to fix the hyper parameters of the network and whenever we simulate, it will first generate a new network and then simulate the covariates.

However, we want to be able to simulate from a fixed network. I.E at most simulate our network according to hyper parameters once, then fix such network, and simulate many different realizations of the covariates and outcomes.

So, I add such functionality by changing the simFromDAG function.

[yiberkeley]

Below is the changed simulation.r:

FormAttributes <- function(DAG, parents, dataform, LTCF, tvals, netind_cl) {
  list(DAG = DAG,
      parents = parents,
      dataform = dataform,
      LTCF = LTCF,
      tvals = tvals,
      netind_cl = netind_cl,
      node_nms = sapply(N(DAG), '[[', "name"),	# node names from the DAG name attributes
      actnodes = attr(DAG, "actnodes"),	# save the action nodes (if defined)
      acttimes = attr(DAG, "acttimes"),	# action time points (if defined)
      attnames = attr(DAG, "attnames"),	# save the attribute names (if defined)
      attrs = attr(DAG, "attrs")			# save the attribute values (if defined)
      )
}
# copy all DAG related attributes from a full data.frame attribute list (skipping row names and such)
# can modify some of the attributes, s.a., dataform, tvals & LTCF if provided
CopyAttributes <- function(attrslist, dataform=NULL, tvals=NULL, LTCF=NULL) {
  list(DAG=attrslist[["DAG"]],
      parents=attrslist[["parents"]],
      dataform = if (!is.null(dataform)) dataform else attrslist[["dataform"]],
      LTCF = if (!is.null(LTCF)) LTCF else attrslist[["LTCF"]],
      tvals = if (!is.null(tvals)) tvals else attrslist[["tvals"]],
      node_nms = attrslist[["node_nms"]],
      actnodes = attrslist[["actnodes"]],
      acttimes = attrslist[["acttimes"]],
      attnames = attrslist[["attnames"]],
      attrs = attrslist[["attrs"]],
      target = attrslist[["target"]]
      )
}

get_allts <- function(obs.df) { # calculate actual t values used in the simulated data (wide format)
  node_nms <- attributes(obs.df)$names[-1]
  all_ts <- as.integer(unique(sapply(strsplit(node_nms, "_"), '[', 2)))
  all_ts <- sort(all_ts[!is.na(all_ts)]) #remove NA/NULL cases and sort
}

get_gennames <- function(obs.df) { # grab generic names (without time points) from simulated data (wide format)
  node_nms <- attributes(obs.df)$names[-1]
  gennames <- as.character(unique(sapply(strsplit(node_nms, "_"), '[', 1)))
  # gennames <- sort(gennames[!is.na(gennames)]) #remove NA/NULL cases and sort
}

# given names of actions, extract the corresponding DAG.action objects from DAG
getactions <- function(DAG, actions) {
  if (!is.character(actions)) stop("actions must be a character vector of action names")
  # check the actions actually exist in the DAG, then grab them
  all_action_DAGs <- A(DAG)
  if (length(all_action_DAGs)<1) stop("DAG objects contains no actions, define actions with DAG+action before simulating action-specific data...")
  all_action_names <- sapply(all_action_DAGs, function(DAG) attr(DAG, "actname"))
  # verify that each action name in actions has only 1 unique match
  action_idx_sel <- NULL
  for (action in actions) {
    action_idx <- which(all_action_names %in% action)
    if (length(action_idx)<1) {
      stop("Couldn't locate action: "%+% action %+% " , first define action by adding it to the DAG object with DAG+action")
    } else if (length(action_idx)>1) {
      stop("Found more than one action matching action name: " %+% action)
    }
    action_idx_sel <- c(action_idx_sel, action_idx)
  }
  actions <- all_action_DAGs[action_idx_sel]
}

# Internal function for simulation data from any DAG
# If prev.data is not NULL all the nodes in DAG will be evaluated in the environment of prev.data alone
simFromDAG <- function(DAG, Nsamp, wide = TRUE, LTCF = NULL, rndseed = NULL, rndseed.reset.node = NULL, prev.data = NULL, verbose = getOption("simcausal.verbose"),oldNet = NULL,oldDat=NULL) {
  assertthat::assert_that(assertthat::is.count(Nsamp) || as.integer(Nsamp)==0L)
  if (!is.null(rndseed)) {
    set.seed(rndseed)
  }
  if (!is.DAG(DAG)) stop("DAG argument must be an object of class DAG")

  # SAMPLING FROM ARBITRARY NODE DISTR FUN:
  sampleNodeDistr <- function(newNodeParams, distr, EFUP.prev, cur.node, expr_str, asis.samp = FALSE) {
    N_notNA_samp <- sum(!EFUP.prev)
    user.env <- cur.node$node.env
    asis.flags <- attr(newNodeParams$dist_params, "asis.flags")
    # print("asis.flags: "); print(asis.flags)
    # if (all(asis.flags)) asis.samp <- TRUE

    standardize_param <- function(distparam) {
      check_len <- function(param) {
        len <- length(param)
        if (len==1) {
          param <- rep.int(param, N_notNA_samp)
        } else if (len==Nsamp) {
          param <- param[!EFUP.prev]
        } else {
          # stop("error while evaluating node "%+% cur.node$name %+%" expression(s): "%+%expr_str%+%".\n One of the distribution parameters evaluated to an incorrect vector length, check syntax.")
          if (verbose) message("evaluating node "%+% cur.node$name %+%" expression(s): "%+%expr_str%+%
              ".\n One of the distribution parameters evaluated to non-standard vector length (its neither 1 nor n), make sure the distribution function knows how to handle it.")
        }
        param
      }

      # print("distparam: "); print(distparam)
      # print("class(distparam): "); print(class(distparam))

      if (is(distparam,"list")) { #if (class(distparam)%in%"list") {
        distparam <- lapply(distparam, check_len)
        distparam <- do.call("cbind", distparam)
      } else if (is.vector(distparam)) {
        distparam <- check_len(distparam)
      } else if (is.matrix(distparam)) {
        if (nrow(distparam)==1) { # for mtx with 1 row -> repeat N_notNA_samp times:
          # distparam <- matrix(distparam, nrow = N_notNA_samp, ncol = ncol(distparam), byrow = TRUE)
        } else if (nrow(distparam)==Nsamp) { # for mtx with Nsamp rows -> subset mtx[!EFUP.prev,]
          distparam <- distparam[!EFUP.prev, , drop=FALSE]
        } else {
          # stop("error while evaluating node "%+% cur.node$name %+%" expression(s): "%+%expr_str%+%".\n One of the distribution parameters evaluated to an incorrect vector length, check syntax.")
          if (verbose) message("evaluating node "%+% cur.node$name %+%
              " expression(s): One of the distribution parameters evaluated to a matrix of non-standard dimensions (its neither 1 nor n rows), make sure the distribution function knows how to handle it.")
        }
      } else {
        # print("expr_str"); print(expr_str)
        # print(class(expr_str))
        # stop("...one of the distribution parameters evaluated to unsported data type...")
        # stop("error while evaluating node "%+% cur.node$name %+%" expression(s): "%+%expr_str%+%".\n One of the formulas evaluated to an unsported data type, check syntax.")
        if (verbose) message("The node " %+%  cur.node$name %+% " expression(s): "%+%expr_str%+%
            ".\n One of the formulas evaluated to unknown data type. Make sure the distribution function knowns how to handle it...")
      }
      distparam
    }

    # 1) Check & get the distribution function from the user-calling environment
    # 2) If it doesn't exist, search the package namespace, if not found throw an exception
    # More efficient version:
    if (is.null(distr.fun <- get0(distr, mode = "function"))) {
      if (!is.null(distr.fun <- get0(distr, envir = user.env, mode = "function"))) {
        if (verbose) message(distr %+% ": note this distribution could not be located in package namespace, simulating from user-defined distribution found under the same name")
      } else {
        namespace_distr <- unlist(strsplit(distr, "::"))
        distr.fun <- tryCatch(getFromNamespace(x = namespace_distr[2], ns = namespace_distr[1], envir = user.env))
        if (inherits(distr.fun, "try-error")) stop(distr %+% ": this node distribution function could not be located")
        # stop(distr %+%": this distribution function can\'t be found")
      }
    }

    # *) Go over distribution parameters and for each parameter:
    # check for consistency; make it an appropriate length; turn each list parameter into matrix;
    dprint("before standardize newNodeParams$dist_params:"); dprint(newNodeParams$dist_params)
    if (!asis.samp) {
      for (idx in seq_along(newNodeParams$dist_params)) {
        if (!asis.flags[[idx]] || (is.vector(newNodeParams$dist_params[[idx]]) && length(newNodeParams$dist_params[[idx]]) == Nsamp))
          newNodeParams$dist_params[[idx]] <- standardize_param(newNodeParams$dist_params[[idx]])
      }
    }
    dprint("after standardize newNodeParams$dist_params:"); dprint(newNodeParams$dist_params)

    # *) Add n argument to dist_params
    newNodeParams$dist_params <-  append(list(n=N_notNA_samp), newNodeParams$dist_params)

    # *) Call the distribution function, passing the distribution parameters:
    newVar_Samp <- try(do.call(distr.fun, newNodeParams$dist_params))

    if(inherits(newVar_Samp, "try-error")) {
      stop("failed while sampling node '" %+% cur.node$name %+% "', from distribution '" %+% distr %+% "'")
    }

    if (asis.samp) {
      return(newVar_Samp)
    # output is a vector
    } else if (is.null(dim(newVar_Samp))||is.vector(newVar_Samp)) {
      newVar <- rep.int(NA, Nsamp)
      if (is.factor(newVar_Samp)) {
        newVar[!EFUP.prev] <- as.numeric(levels(newVar_Samp))[newVar_Samp]
        newVar <- factor(newVar)
      } else {
        newVar[!EFUP.prev] <- newVar_Samp
      }
      return(newVar)
    # output is a matrix -> test dimensionality is correct (equal length(names)), but do nothing else
    } else if (is.matrix(newVar_Samp)) {
      if (ncol(newVar_Samp)!=length(cur.node$mv.names)) stop("dimensionality (number of columns) sampled by multivariate node is not equal to length of `name` argument for node: " %+% cur.node$name)
      colnames(newVar_Samp) <- cur.node$mv.names
      return(newVar_Samp)
    } else {
      stop("unrecognized sampled variable output type, only vectors or matrices are allowed, node: " %+% cur.node$name)
    }
  }

  #---------------------------------------------------------------------------------
  # Iteratively create observed DAG nodes (wide format data.fame of Nsamp observations)
  #---------------------------------------------------------------------------------
  EFUP.prev <- rep(FALSE, Nsamp)
  LTCF.prev <- rep(FALSE, Nsamp)

  obs.df <- data.frame(ID = seq(vector(length = Nsamp))) # obs.df <- data.frame(ID = seq(1:Nsamp))
  obs.dt <- data.table(obs.df)

  # for efficiency, allocate columns in DT to maximum number of variables, if it exceeds 1025 (default)
  DT.cols <- (length(DAG)+1)
  if (DT.cols > 1025) {
    alloc.col(obs.dt, DT.cols)
  }
  # alloc.col(obs.dt, max(min(length(DAG)+1, 1000), 200))

  #---------------------------------------------------------------------------------
  # CHECKS PERFORMED DURING EVALUTION:
  #---------------------------------------------------------------------------------
  # at each node check:
  # *) (formulas) check all formulas are evaluable R expressions (try(eval))
  # *) (formulas) check each cur.node formula references:
  #	node names that already exist, and all node names must have order < cur.node$order

  #---------------------------------------------------------------------------------
  # consider pre-allocating memory for newNodeParams and newVar
  #---------------------------------------------------------------------------------
  newVar <- vector(length = Nsamp)
  NvarsPerNode <- sapply(DAG, function(node) length(node[["mv.names"]]))
  totalNvars <- sum(NvarsPerNode)
  MVvars <- unlist(lapply(DAG, function(node) node[["mv.names"]]))
  MVvarMapNode <- NULL
  for (node in DAG) {
    for (mv.name in node[["mv.names"]]) MVvarMapNode <- append(MVvarMapNode, list(list(MVname = mv.name, Node = node$name)))
  }
  names(MVvarMapNode) <- MVvars

  NodeParentsNms <- vector("list", length = length(DAG))	# list that will have parents names for each node
  names(NodeParentsNms) <- names(DAG)

  # initialize the node formula evaluator class (Define_sVar):
  node_evaluator <- Define_sVar$new() # netind_cl = netind_cl
  t_pts <- NULL
  netind_cl <- NULL

  for (cur.node in DAG) {
    t <- cur.node$t # define variable for a current time point t
    t_new <- !(t %in% t_pts) # set to TRUE when switch to new time point t occurs otherwise FALSE
    t_pts <- c(t_pts, t)
    t_pts <- as.integer(unique(t_pts)) # vector that keeps track of unique timepoints t
    gnodename <- as.character(unlist(strsplit(cur.node$name, "_"))[1])

    #------------------------------------------------------------------------
    # CHECKS NODE DISTRIBUTIONS & EVALUATE NODE DISTRIBUTION PARAMS PROVIDED IN NODE FORMULA(S)
    #------------------------------------------------------------------------
    # TO DO:
    # Need to make sure that action attributes never get separate columns in df, since they are just constants
    # The attributes should be added to the eval env as variables
    # See if current expr naming strucutre in Define_sVar can be recycled for sampling from multivariate densities
    # NEED TO ADDRESS: WHEN A CONSTANT (length==1) IS PASSED AS A PARAMETER, DON'T ALWAYS WANT TO TURN IT INTO A VECTOR OF LENGTH n
    #------------------------------------------------------------------------
    # setting the node-specific user calling environment for the evaluator:
    # node_evaluator$set.user.env(cur.node$node.env)
    eval_expr_res <- node_evaluator$eval.nodeforms(cur.node = cur.node, data.df = if (!is.null(prev.data)) prev.data else obs.dt)
    par.names <- unique(unlist(lapply(eval_expr_res, '[[', "par.nodes")))
    eval_dist_params <- lapply(eval_expr_res, '[[' ,"evaled_expr")
    attr(eval_dist_params, "asis.flags") <- attr(cur.node$dist_params, "asis.flags")
    newNodeParams <- list(dist_params = eval_dist_params, par.names = par.names)

    # dprint("cur.node t:" %+% cur.node$t)
    # dprint("original nodeform expr:"); dprint(lapply(cur.node$dist_params, function(expr) as.character(expr)))
    # dprint("full eval_expr_res:"); dprint(eval_expr_res)
    # dprint("eval_dist_params:"); dprint(eval_dist_params)

    #------------------------------------------------------------------------
    # SAMPLE NEW COVARIATE BASED ON DISTR PARAMS
    #------------------------------------------------------------------------

    # reset the seed if the current node matches argument rndseed.reset.node
    if (!is.null(rndseed) && !is.null(rndseed.reset.node)) {
      if (rndseed.reset.node %in% cur.node$name) set.seed(NULL)
    }

  
    distr <- cur.node$distr
    
    if(is.null(oldNet)){
      
      
    if ("DAG.net" %in% class(cur.node)) {

      if (!is.null(netind_cl) && verbose) message("Network is being re-sampled during the same simulation!")
      message("1,Network is being re-sampled ")
      distr <- cur.node$netfun

      # ***************************************************************
      # USING AS-IS FOR SAMPLING NETWORK BECAUSE THE OUTPUT IS A MATRIX. THIS IS A BAD WAY TO SOLVE THIS.
      # A BETTER SOLUTION IS TO ALLOW THE distr RESULT TO BE A VECTOR OR MATRIX (for multivariate RVs)
      # *****************************************************************
      NetIndMat <- sampleNodeDistr(newNodeParams = newNodeParams, distr = distr, EFUP.prev = EFUP.prev,
                                  cur.node = cur.node, expr_str = cur.node$dist_params, asis.samp = TRUE)
      cur.node$Kmax <- ncol(NetIndMat)
      netind_cl <- NetIndClass$new(nobs = Nsamp, Kmax = cur.node$Kmax)
      netind_cl$NetInd <- NetIndMat
      node_evaluator$set.net(netind_cl) # set the network for the evaluator:

      # Kmax.new <- ncol(NetIndMat)
      # if (as.integer(Kmax.new) != as.integer(cur.node$Kmax)) {
      #   if (verbose) message("Kmax is reset to a new value; old Kmax: " %+% cur.node$Kmax %+% "; new Kmax: " %+% Kmax.new)
      #   cur.node$Kmax <- Kmax.new
      # }
      # node_evaluator$Kmax <- netind_cl$Kmax
      newVar <- NULL
      NodeParentsNms <- NodeParentsNms[-which(names(NodeParentsNms) %in% cur.node$name)]

    } else {
      # print("cur.node: "); print(cur.node$name)
      # print("newNodeParams$par.names: "); print(newNodeParams$par.names)
      # print("MVvarMapNode"); print(MVvarMapNode)
      if (!is.null(newNodeParams$par.names)) {
        for (idx in seq_along(newNodeParams$par.names)) {
          choose_node_idx <- names(MVvarMapNode) %in% newNodeParams$par.names[idx]
          if (sum(choose_node_idx)>0) newNodeParams$par.names[idx] <- MVvarMapNode[choose_node_idx][[1]][["Node"]]
        }
        NodeParentsNms[[cur.node$name]] <- c(NodeParentsNms[[cur.node$name]], newNodeParams$par.names)
      }

      newVar <- sampleNodeDistr(newNodeParams = newNodeParams, distr = distr, EFUP.prev = EFUP.prev,
                                cur.node = cur.node, expr_str = cur.node$dist_params)
    }
      
      
      
      
      
    }
    else
    {
      
      
      
      
      if ("DAG.net" %in% class(cur.node)) {
        
      
        message("Using the old network structure, but just generating covariate")
        distr <- cur.node$netfun
        
        NetIndMat <- attributes(oldDat)$netind_cl$NetInd
        cur.node$Kmax <- attributes(oldDat)$netind_cl$Kmax
        netind_cl <- attributes(oldDat)$netind_cl
        netind_cl$NetInd <- NetIndMat
        node_evaluator$set.net(netind_cl) # set the network for the evaluator:
      
        
        newVar <- NULL
        NodeParentsNms <- NodeParentsNms[-which(names(NodeParentsNms) %in% cur.node$name)]
        
      } else {
  
        if (!is.null(newNodeParams$par.names)) {
          for (idx in seq_along(newNodeParams$par.names)) {
            choose_node_idx <- names(MVvarMapNode) %in% newNodeParams$par.names[idx]
            if (sum(choose_node_idx)>0) newNodeParams$par.names[idx] <- MVvarMapNode[choose_node_idx][[1]][["Node"]]
          }
          NodeParentsNms[[cur.node$name]] <- c(NodeParentsNms[[cur.node$name]], newNodeParams$par.names)
        }
        
        newVar <- sampleNodeDistr(newNodeParams = newNodeParams, distr = distr, EFUP.prev = EFUP.prev,
                                  cur.node = cur.node, expr_str = cur.node$dist_params)
      }
      
    }
    

    if (!is.null(newVar)) {
      # already performed inside sampleNodeDistr, so commented out:
      # newVar[EFUP.prev] <- NA
      # uncomment to assign 0 instead of NA for missing (after EOF=TRUE node evals to 1)
      # newVar[EFUP.prev] <- 0
      #------------------------------------------------------------------------
      # LTCF is a generic name for a failure/censoring node, after failure occurred (node=1), carry forward (impute) all future covariate values (for example, survival outcome with EFUP=TRUE)
      #------------------------------------------------------------------------
      if (!is.null(t)) {	# check time points were specified by user
        # check current t is past the first time-point & we want to last time-point carried forward (LTCF)
        if ((sum(LTCF.prev) > 0) & (t > t_pts[1]) & !is.null(LTCF)) {
          # Carry last observation forward for those who reached EFUP (failed or but not censored)
          prevtVarnm <- gnodename %+% "_" %+% (t-1)
          # Check first that the variable exists (has been defined at the previous time point), if it doesn't exist, just assign NA
          # if(!any(names(obs.df)%in%prevtVarnm)) {
          if(!any(names(obs.dt)%in%prevtVarnm)) {
            warning(gnodename%+%": is undefined at t="%+%(t-1)%+%", hence cannot be carried forward to t="%+%t)
            newVar[LTCF.prev] <- NA
          } else {
            # newVar[LTCF.prev] <- obs.df[LTCF.prev, (gnodename %+% "_" %+% (t-1))]
            # newVar[LTCF.prev] <- obs.dt[LTCF.prev, (gnodename %+% "_" %+% (t-1)), with = FALSE]
            newVar[LTCF.prev] <- obs.dt[LTCF.prev, ][[(gnodename %+% "_" %+% (t-1))]]
          }
        }
        # carry forward the censoring indicator as well (after evaluating to censored=1) - disabled call doLTCF() with censoring node name instead
        # if ((!is.null(LTCF))) {
        # 	if (is.EFUP(cur.node) & (!(LTCF%in%gnodename))) {
        # 		newVar[EFUP.prev & !(LTCF.prev)] <- 1
        # 	}
        # }
      }
      #------------------------------------------------------------------------
      # modify the observed data by adding new sampled covariate
      # *** TODO: need to internally save the evaluation results EFU.flag for the following lines to work correctly:
        # in eval.target():
        # vec_EFUP <- sapply(N(DAG)[outnode_nms], is.EFUP)
        # in eval.MSM():
        # vec_EFUP <- sapply(DAG[outnode_nms], is.EFUP)
        # in doLTCF():
        # if  (is.EFUP(cur.node)&(cur.outcome))
      #------------------------------------------------------------------------

      # 2 ways to assign several new vars with data.table (1st is the fastest)
      if (is.matrix(newVar)) {
        for (col in cur.node$mv.names) {
          obs.dt[!EFUP.prev, (col) := newVar[, col]]
          # obs.df <- within(obs.df, {assign(col, newVar[, col])})
        }
      } else {
        obs.dt[, (cur.node$name) := newVar]
        # obs.df <- within(obs.df, {assign(cur.node$name, newVar)})
      }
      # another version, requires making a copy of newVar
      # obs.dt[,(cur.node$name):=as.data.table(newVar)]

      if (!is.null(cur.node$EFU)) {
        # EFU is now evaluated separately for each observation, only when evaluates to TRUE then the current node starts acting as a censoring variable
        # EFU.flag <- node_evaluator$eval.EFU(cur.node = cur.node, data.df = if (!is.null(prev.data)) prev.data else obs.df)
        EFU.flag <- node_evaluator$eval.EFU(cur.node = cur.node, data.df = if (!is.null(prev.data)) prev.data else obs.dt)
        EFU.TRUE <- EFU.flag %in% TRUE
        # EFU.TRUE <- which(EFU.flag %in% TRUE)
        # print("EFU.flag: "); print(EFU.flag); print(EFU.flag %in% TRUE)
        # check for censoring only when at least one observation had EFU=TRUE
        if (sum(EFU.TRUE)>0) {
          # EFUP.now <- (obs.df[,ncol(obs.df)]%in%1) & (EFU.TRUE)
          EFUP.now <- (obs.dt[[cur.node$name[1]]]%in%1) & (EFU.TRUE)
          EFUP.prev <- (EFUP.prev | EFUP.now)
          if ((!is.null(LTCF)) && (LTCF%in%gnodename[1]) && is.null(prev.data)) { # is this node the one to be carried forward (LTCF node)? mark the observations with value = 1 (carry forward)
            LTCF.prev <- (LTCF.prev | EFUP.now)
          }
        }
      }
    }
  }

  # Collect all attributes to be assigned to the obs data
  # all_ts <- get_allts(obs.df) # calculate actual time values used in the simulated data
  all_ts <- get_allts(obs.dt) # calculate actual time values used in the simulated data

  if (sum(LTCF.prev)>0) {	# only change the LTCF attribute if outcome carry forward imputation was really carried out for at least one obs
    LTCF_flag <- LTCF
  } else {
    LTCF_flag <- NULL
  }

  # -----------------------------------------------------------------------------------
  # Remove latent variables from the simulated dataset (if any were declared as latent)
  # -----------------------------------------------------------------------------------
  latent.v <- attr(DAG, "latent.v")
  if (!is.null(latent.v)) {
    # excl.cols <- colnames(obs.df)%in%latent.v
    excl.cols <- colnames(obs.dt)%in%latent.v
    if (sum(excl.cols) > 0) {
      # obs.df <- obs.df[,!excl.cols]
      obs.dt[,(colnames(obs.dt)[excl.cols]) := NULL]
    }
  }

  newattrs <- FormAttributes(DAG = DAG, parents = NodeParentsNms, dataform = "wide", LTCF = LTCF_flag, tvals = all_ts, netind_cl = netind_cl)
  newattrs_names <- names(newattrs)
  for (attr_name in newattrs_names) {
    setattr(x = obs.dt, name = attr_name, value = newattrs[[attr_name]])
  }

  if (!wide) {
    if (length(all_ts)>1) {		# only perform conversion when there is more than one t in the simulated data
      # obs.df <- DF.to.long(obs.df)
      obs.dt <- DF.to.longDT(obs.dt)
    } else {
      warning("Simulated data returned in wide format. Can't convert to long format, since only one time-point was detected.")
    }
  }
  class(obs.dt) <- "data.frame"
  # reset the seed to NULL if it was previously set
  if (!is.null(rndseed)) {
    set.seed(NULL)
  }
  # return(obs.df)
  return(obs.dt)
}

#' Simulate Observed Data
#'
#' This function simulates observed data from a DAG object.
#' @param DAG A DAG objects that has been locked with set.DAG(DAG). Observed data from this DAG will be simulated.
#' @param n Number of observations to sample.
#' @param wide A logical, if TRUE the output data is generated in wide format, if FALSE, the output longitudinal data in generated in long format
#' @param LTCF If forward imputation is desired for the missing variable values, this argument should be set to the name of the node that indicates the end of follow-up event. See the vignette, \code{\link{sim}} and \code{\link{doLTCF}} for additional details.
#' @param rndseed Seed for the random number generator.
#' @param rndseed.reset.node When \code{rndseed} is specified, use this argument to specify the name of the \code{DAG} node at which the random number generator seed is reset back to \code{NULL} (simulation function will call \code{set.seed(NULL)}).
#' Can be useful if one wishes to simulate data using the set seed \code{rndseed} only for the first K nodes of the DAG and use an entirely random sample when simulating the rest of the nodes starting at K+1 and on.
#' The name of such (K+1)th order \code{DAG} node should be then specified with this argument.
#' @param verbose Set to \code{TRUE} to print messages on status and information to the console.
#'  Turn this off by default using options(simcausal.verbose=FALSE).
#' @return A \code{data.frame} where each column is sampled from the conditional distribution specified by the corresponding \code{DAG} object node.
#' @family simulation functions
#' @seealso \code{\link{simfull}} - a wrapper function for simulating full data only; \code{\link{sim}} - a wrapper function for simulating both types of data; \code{\link{doLTCF}} for forward imputation of the missing values in already simulating data; \code{\link{DF.to.long}}, \code{\link{DF.to.longDT}} - converting longitudinal data from wide to long formats.
#' @export
simobs <- function(DAG, n, wide = TRUE, LTCF = NULL, rndseed = NULL, rndseed.reset.node = NULL, verbose = getOption("simcausal.verbose"),oldNet = NULL,oldDat=NULL) {
  if (!is.DAG(DAG)) stop("DAG argument must be an object of class DAG")
  simFromDAG(DAG=DAG, Nsamp=n, wide=wide, LTCF=LTCF, rndseed=rndseed, rndseed.reset.node=rndseed.reset.node, verbose=verbose,oldNet = oldNet,oldDat=oldDat)
}

#' Simulate Full Data (From Action DAG(s))
#'
#' This function simulates full data based on a list of intervention DAGs, returning a list of \code{data.frame}s.
#' @param actions Actions specifying the counterfactual DAG. This argument must be either an object of class DAG.action or a list of DAG.action objects.
#' @param n Number of observations to sample.
#' @param wide A logical, if TRUE the output data is generated in wide format, if FALSE, the output longitudinal data in generated in long format
#' @param LTCF If forward imputation is desired for the missing variable values, this argument should be set to the name of the node that indicates the end of follow-up event. See the vignette, \code{\link{sim}} and \code{\link{doLTCF}} for additional details.
#' @param rndseed Seed for the random number generator.
#' @param rndseed.reset.node When \code{rndseed} is specified, use this argument to specify the name of the \code{DAG} node at which the random number generator seed is reset back to \code{NULL} (simulation function will call \code{set.seed(NULL)}).
#' Can be useful if one wishes to simulate data using the set seed \code{rndseed} only for the first K nodes of the DAG and use an entirely random sample when simulating the rest of the nodes starting at K+1 and on.
#' The name of such (K+1)th order \code{DAG} node should be then specified with this argument.
#' @param verbose Set to \code{TRUE} to print messages on status and information to the console.
#'  Turn this off by default using options(simcausal.verbose=FALSE).
#' @return A named list, each item is a \code{data.frame} corresponding to an action specified by the actions argument, action names are used for naming these list items.
#' @family simulation functions
#' @seealso \code{\link{simobs}} - a wrapper function for simulating observed data only; \code{\link{sim}} - a wrapper function for simulating both types of data; \code{\link{doLTCF}} for forward imputation of the missing values in already simulating data; \code{\link{DF.to.long}}, \code{\link{DF.to.longDT}} - converting longitudinal data from wide to long formats.
#' @export
simfull <- function(actions, n, wide = TRUE, LTCF = NULL, rndseed = NULL, rndseed.reset.node = NULL, verbose = getOption("simcausal.verbose"),oldNet = NULL,oldDat=NULL) {
  if (!is.null(rndseed)) {
    set.seed(rndseed)
  }
  if (("DAG.action" %in% class(actions))) {
    actname <- get.actname(actions)
    actions <- list(actions)
    names(actions) <- actname
  } else if ("list" %in% class(actions)) {
    checkDAGs <- sapply(actions, is.DAG.action)
    if (!all(checkDAGs)) stop("argument 'actions': all elements of the list must be DAG.action objects")
  } else {
    stop("argument 'actions' must be an object of class 'DAG.action' or a list of 'DAG.action' objects")
  }

  fulldf.list <- list()	# list that will contain full data (each item = one dataframe per intervention)
  for (idx_act in c(1:length(actions))) {	# loop over each action/intervention/regimen DAG
    actname <- names(actions)[idx_act]
    actDAG <- actions[[idx_act]]
    fulldf <- try(simFromDAG(DAG=actDAG, Nsamp=n, wide=wide, LTCF=LTCF, rndseed=rndseed, rndseed.reset.node=rndseed.reset.node, verbose=verbose,oldNet = oldNet ,oldDat=oldDat))
    if(inherits(fulldf, "try-error")) {
    	stop("simulating full data for action '"%+%actname%+%"' resulted in error...", call. = FALSE)
    }
    # adding action indicator column to the data.frame
    # newattrs <- CopyAttributes(attrslist=attributes(fulldf)) 	# save attributes
    # actindvar <- rep.int(actname, times=nrow(fulldf))
    # fulldf <- cbind(action=actindvar, fulldf)
    # attributes(fulldf) <- c(attributes(fulldf), newattrs)
    # attr(fulldf, "node_nms") <- c(actname, attr(fulldf, "node_nms"))
    # names(attr(fulldf, "node_nms"))[1] <- actname
    fulldf.list <- c(fulldf.list, list(fulldf))
  }
  if (length(actions)!=length(fulldf.list)) stop("simfull(): couldn't create separate df for each action DAG")
  names(fulldf.list) <- names(actions)
  return(fulldf.list)
}

#' Simulate Observed or Full Data from \code{DAG} Object
#'
#' This function simulates full data based on a list of intervention DAGs, returning a list of \code{data.frame}s. See the vignette for examples and detailed description.
#' @section Forward Imputation:
#' By default, when LTCF is left unspecified, all variables that follow after any end of follow-up (EFU) event are set to missing (NA).
#' The end of follow-up event occurs when a binary node of type \code{EFU=TRUE} is equal to 1, indicating a failing or right-censoring event.
#' To forward impute the values of the time-varying nodes after the occurrence of the \code{EFU} event, set the LTCF argument to a name of the EFU node representing this event.
#' For additional details and examples see the vignette and \code{\link{doLTCF}} function.
#' @param DAG A DAG objects that has been locked with set.DAG(DAG). Observed data from this DAG will be simulated if actions argument is omitted.
#' @param actions Character vector of action names which will be extracted from the DAG object. Alternatively, this can be a list of action DAGs selected with \code{A(DAG)} function, in which case the argument \code{DAG} is unused. When \code{actions} is omitted, the function returns simulated observed data (see \code{simobs}).
#' @param n Number of observations to sample.
#' @param wide A logical, if TRUE the output data is generated in wide format, if FALSE, the output longitudinal data in generated in long format
#' @param LTCF If forward imputation is desired for the missing variable values, this argument should be set to the name of the node that indicates the end of follow-up event.
#' @param rndseed Seed for the random number generator.
#' @param rndseed.reset.node When \code{rndseed} is specified, use this argument to specify the name of the \code{DAG} node at which the random number generator seed is reset back to \code{NULL} (simulation function will call \code{set.seed(NULL)}).
#' Can be useful if one wishes to simulate data using the set seed \code{rndseed} only for the first K nodes of the DAG and use an entirely random sample when simulating the rest of the nodes starting at K+1 and on.
#' The name of such (K+1)th order \code{DAG} node should be then specified with this argument.
#' @param verbose Set to \code{TRUE} to print messages on status and information to the console.
#'  Turn this off by default using options(simcausal.verbose=FALSE).
#' @return If actions argument is missing a simulated data.frame is returned, otherwise the function returns a named list of action-specific simulated data.frames with action names giving names to corresponding list items.
#' @family simulation functions
#' @seealso \code{\link{simobs}} - a wrapper function for simulating observed data only; \code{\link{simfull}} - a wrapper function for simulating full data only; \code{\link{doLTCF}} - forward imputation of the missing values in already simulating data; \code{\link{DF.to.long}}, \code{\link{DF.to.longDT}} - converting longitudinal data from wide to long formats.
#' @example tests/examples/sim.impute.examples12.R
#' @references Sofrygin O, van der Laan MJ, Neugebauer R (2017).
#' "simcausal R Package: Conducting Transparent and Reproducible Simulation Studies of Causal Effect Estimation with Complex Longitudinal Data."
#' Journal of Statistical Software, 81(2), 1-47. doi: 10.18637/jss.v081.i02.
#' @export
sim <- function(DAG, actions, n, wide = TRUE, LTCF = NULL, rndseed = NULL, rndseed.reset.node = NULL, verbose = getOption("simcausal.verbose"),oldNet = NULL,oldDat=NULL) {
  # *) check if actions argument is missing -> simulate observed data from the DAG
  # *) if actions consist of characters, try to extract those actions from the DAG and simulate full data
  if (!is.integerish(n)) stop("n must be an integer")
  # if (!(n >= 1)) stop("n must be a positive integer")

  # SIMULATE OBSERVED DATA FROM DAG (if no actions)
  if (missing(actions)) {
    if (!is.DAG(DAG)) stop("DAG argument must be an object of class DAG")
    if (!is.DAGlocked(DAG)) stop("call set.DAG() before attempting to simulate data from DAG")
    if (verbose) message("simulating observed dataset from the DAG object")
    return(simobs(DAG = DAG, n = n, wide = wide, LTCF = LTCF, rndseed = rndseed, rndseed.reset.node = rndseed.reset.node, verbose = verbose,oldNet = oldNet,oldDat=oldDat))
  # SIMULATE FULL DATA FROM ACTION DAGs (when actions are present)
  } else {
    if (is.character(actions)) { # grab appropriate actions from the DAG by name
      if (!is.DAG(DAG)) stop("DAG argument must be an object of class DAG")
      actions <- getactions(DAG, actions)
    } else if (is.list(actions)) {
      checkDAGs <- sapply(actions, is.DAG)
      if (!all(checkDAGs)) stop("actions must be either a list of DAGs or a vector of action names")
    } else {
      stop("argument actions must be a character vector of action names or a list of action DAGs")
    }
    if (verbose) message("simulating action-specific datasets for action(s): "%+%paste(names(actions), collapse = " "))
    return(simfull(actions = actions, n = n, wide = wide, LTCF = LTCF, rndseed = rndseed, rndseed.reset.node = rndseed.reset.node, verbose = verbose,oldNet = oldNet,oldDat=oldDat))
  }
}

#' Missing Variable Imputation with Last Time Point Value Carried Forward (LTCF)
#'
#' Forward imputation for missing variable values in simulated data after a particular end of the follow-up event. The end of follow-up event is defined by the node of type \code{EOF=TRUE} being equal to 1.
#' @section Details:
#' The default behavior of the \code{sim} function consists in setting all nodes that temporally follow an \code{EFU} node whose simulated value is 1 to missing (i.e., \code{NA}).
#' The argument \code{LTCF} of the \code{sim} function can however be used to change this default behavior and impute some of these missing values with \emph{last time point value carried forward} (LTCF).
#' More specifically, only the missing values of time-varying nodes (i.e., those with non-missing \code{t} argument) that follow the end of follow-up event encoded by the \code{EFU} node specified by the \code{LTCF} argument will be imputed.
#' One can use the function \code{doLTCF} to apply the \emph{last time point value carried forward} (LTCF) imputation to an existing simulated dataset obtained from the function \code{sim} that was called with its default imputation setting (i.e., with no \code{LTCF} argument).
#' Illustration of the use of the LTCF imputation functionality are provided in the package vignette.
#'
#' The first example below shows the default data format of the \code{sim} function after an end of the follow-up event and how this behavior can be modified to generate data with LTCF imputation by either using the \code{LTCF} argument of the
#' \code{sim} function or by calling the \code{doLTCF} function. The second example demonstrates how to use the \code{doLTCF} function to perform LTCF imputation on already existing data simulated with the \code{sim} function based on its default non-imputation behavior.
#' @param data Simulated \code{data.frame} in wide format
#' @param LTCF Character string specifying the outcome node that is the indicator of the end of follow-up (observations with value of the outcome variable being 1 indicate that the end of follow-up has been reached). The outcome variable must be a binary node that was declared with \code{EFU=TRUE}.
#' @return Modified \code{data.frame}, all time-varying missing variables after the \code{EFU} outcome specified in \code{LTCF} are forward imputed with their last available non-missing value.
#' @family data manipulation functions
#' @seealso \code{\link{sim}}, \code{\link{simobs}} and \code{\link{simfull}} for simulating data with and without carry forward imputation.
#' @example tests/examples/sim.impute.examples12.R
#' @export
doLTCF <- function(data, LTCF) {
  DAG <- attr(data, "DAG")
  Nsamp <- nrow(data)
  LTCF.prev <- rep(FALSE,Nsamp)
  t_pts <- NULL
  cur.outcome <- FALSE
  if (is.longfmt(data)) stop("Last Timepoint Carry Forward (LTCF) can only be executed for data in wide format")
  #------------------------------------------------------------------------
  # Loop over nodes in DAG and keep track of EFUP observations in the data (already sampled)
  #------------------------------------------------------------------------
  # all_node_names <- as.character(unlist(lapply(DAG, "[[", "mv.names")))
  # [[1]][["mv.names"]]
  # for (cur.node in DAG) {
  #   print(cur.node[["name"]])
  # }
  for (cur.node in DAG) {
  # for (cur.node %in% all_node_names)
    # gnodename <- as.character(unlist(strsplit(cur.node$name, "_"))[1])
    # cur.node.name <- cur.node[["mv.names"]]
    gnodename <- as.character(unlist(lapply(strsplit(cur.node[["mv.names"]], "_"), "[[", 1)))
    # gnodename <- as.character(unlist(strsplit(cur.node[["mv.names"]], "_"))[1])
    # gnodename <- as.character(unlist(strsplit(cur.node, "_"))[1])

    t <- cur.node$t # define variable for a current time point t
    t_new <- !(t%in%t_pts) # set to TRUE when switch to new time point t occurs otherwise FALSE
    t_pts <- c(t_pts,t); t_pts <- as.integer(unique(t_pts)) # vector that keeps track of unique timepoints t
    if (gnodename[1]%in%LTCF) { # if this generic node name is the outcome we need to carry forward this observation when it evaluates to 1
      cur.outcome <- TRUE
    } else {
      cur.outcome <- FALSE
    }
    if (!is.null(t)) {  # check time points were specified by user
      # check current t is past the first time-point & we want to last time-point carried forward (LTCF)
      if ((sum(LTCF.prev) > 0) & (t > t_pts[1])) { # Carry last observation forward for those who reached EFUP (outcome value = 1)
        if (!gnodename[1]%+%"_"%+%(t-1) %in% names(data)){
          stop("Cannot perform last time point carry forward (LTCF) imputation. No node name '" %+% gnodename%+%"_"%+%(t-1) %+% "' found. Please make sure to define appropriate time-point for this node.")
        }
        for (name in gnodename)
          data[LTCF.prev, name%+%"_"%+%(t)] <- data[LTCF.prev, (name%+%"_"%+%(t-1))]
      }
    }
    if  (is.EFUP(cur.node)&(cur.outcome)) { # if cur.node is EFUP=TRUE type set all observations that had value=1 to LTCF.prev[indx]=TRUE
      LTCF.prev <- (LTCF.prev | (data[,cur.node$name]%in%1))
    }
    # browser()
    # if (!is.null(cur.node[["mv.names"]])) {
    #   for (cur.node in cur.node[["mv.names"]]) {
    #   }
    # }
  }



  if (sum(LTCF.prev)>0) {	# only change the LTCF attribute if outcome carry forward imputation was carried out for at least one obs
    LTCF_flag <- LTCF
  } else {
    LTCF_flag <- NULL
  }
  attr(data, "LTCF") <- LTCF_flag
  # newattrs <- CopyAttributes(attrslist=attrslist, LTCF=LTCF_flag) 	# save attributes from input format data
  # attributes(simdf_long) <- c(attributes(simdf_long), newattrs)
  return(data)
}

#' Convert Data from Wide to Long Format Using \code{reshape}
#'
#' This utility function takes a simulated data.frame in wide format as an input and converts it into a long format (slower compared to \code{\link{DF.to.longDT}}).
#'
#' Keeps all covariates that appear only once and at the first time-point constant (carry-forward).
#'
#' All covariates that appear fewer than range(t) times are imputed with NA for missing time-points.
#'
#' Observations with all NA's for all time-varying covariates are removed.
#'
#' When removing NA's the time-varying covariates that are attributes (attnames) are not considered.
#'
#' @param df_wide A \code{data.frame} in wide format
#' @return A \code{data.frame} object in long format
#' @seealso \code{\link{DF.to.longDT}} - a faster version of \code{DF.to.long} that uses \code{data.table} package
#' @family data manipulation functions
#' @export
DF.to.long <- function(df_wide) {
  if (is.data.table(df_wide)) class(df_wide) <- "data.frame"

  Nsamp <- nrow(df_wide)
  all_ts <- get_allts(df_wide) # calculate actual time values used in the simulated data
  attnames <- attr(df_wide, "attnames")
  node_nms <- attr(df_wide, "node_nms")
  node_nms_split <- strsplit(node_nms, "_")
  node_nms_vec <- sapply(node_nms_split, '[', 1)
  node_nms_unq <- unique(node_nms_vec)

  # if there are no time-points (t) attributes, then the long vs. wide format is undefined.
  if (length(all_ts)==0) return(df_wide)

  varying <- list()
  v.names <- NULL
  bsl.names <- NULL
  for (cur_node_name in node_nms_unq) {
    idx_node <- sapply(node_nms_split, function(t_nodename) t_nodename[1]%in%cur_node_name)
    count_ts <- sum(idx_node)
    node_t_vals <- unlist(sapply(node_nms_split, function(t_nodename) if(t_nodename[1]%in%cur_node_name) as.integer(t_nodename[2])))
    # if node name has no "_", assume its bsl and assign first time-point:
    if (all(is.na(node_t_vals))) node_t_vals <- all_ts[1]

    dprint("cur_node_name: "%+%cur_node_name);
    dprint("all_ts"); dprint(all_ts)
    # dprint("idx_node"); dprint(idx_node)
    dprint("count_ts"); dprint(count_ts)
    dprint("node_t_vals"); dprint(node_t_vals)

    if (count_ts==length(all_ts)) { # the number of time-points is equal to all time-points => # everything is well defined
      varying <- c(varying, list(node_nms[idx_node]))
      v.names <- c(v.names, cur_node_name)
    } else if ((count_ts==1) & (all_ts[all_ts%in%node_t_vals][1]==all_ts[1])) { # only one time-point and its the first-one
      # node is baseline: everything is fine, rename baseline var removing the _t
      bsl.dat.name <- node_nms[which(node_nms_vec%in%cur_node_name)][1]
      bsl.names <- c(bsl.names, cur_node_name)
      if (bsl.dat.name!=cur_node_name) {
        names(df_wide)[names(df_wide) %in% bsl.dat.name] <- cur_node_name
      }
    } else { # in between 1 & 2 => need to impute to length(all_ts) filling up with NAs
      ts_to_NAs <- all_ts[!all_ts%in%node_t_vals]
      ts_not_NAs <- all_ts[all_ts%in%node_t_vals]
      NA_nodenames <- cur_node_name%+%"_"%+%ts_to_NAs # define these node names as missing and add them to df
      df_NAts <- matrix(NA, nrow=Nsamp, ncol=length(ts_to_NAs))
      colnames(df_NAts) <- NA_nodenames
      df_wide <- cbind(df_wide, df_NAts)
      all_vnames <- cur_node_name%+%"_"%+%all_ts
      varying <- c(varying, list(all_vnames))
      v.names <- c(v.names, cur_node_name)
    }
  }
  simdf_long <- reshape(data=df_wide, varying=varying, v.names=v.names, timevar = "t", times = all_ts, sep="_", direction="long")
  simdf_long <- simdf_long[order(simdf_long$ID, simdf_long$t),] # 1) long format is sorted by t first, then ID -> needs to be sorted by ID then time
  simdf_long <- simdf_long[,-ncol(simdf_long)]  # 2) remove last column and reorder columns ID, t, ...
  if (length(attnames)>0) v.names <- v.names[!v.names%in%attnames]
  dprint("v.names"); dprint(v.names)
  # dprint("simdf_long after reshape"); dprint(simdf_long)
  simdf_long <- simdf_long[rowSums(is.na(simdf_long[,v.names]))<length(v.names),] # 3) remove last NA rows (after censoring or EFUP)
  newattrs <- CopyAttributes(attrslist=attributes(df_wide), dataform="long") 	# save attributes from input wide format data
  attributes(simdf_long) <- c(attributes(simdf_long), newattrs)
  attr(simdf_long, "v.names") <- v.names
  simdf_long
}

#' @import data.table
NULL
#' Faster Conversion of Data from Wide to Long Format Using \code{dcast.data.table}
#'
#' Faster utility function for converting wide-format \code{data.frame} into a long format.
#' Internally uses \pkg{data.table} package functions \code{melt.data.table} and \code{dcast.data.table}.
#'
#' Keeps all covariates that appear only once and at the first time-point constant (carry-forward).
#'
#' All covariates that appear fewer than range(t) times are imputed with NA for missing time-points.
#'
#' Observations with all NA's for all time-varying covariates are removed.
#'
#' When removing NA's the time-varying covariates that are attributes (attnames) are not considered.
#'
#' @param df_wide A \code{data.frame} or \code{data.table} in wide format
#' @param return_DF \code{TRUE} (default) to return a \code{data.frame}, \code{FALSE} returns a \code{data.table}
#' @return A \code{data.frame} in long format
#' @family data manipulation functions
#' @export
DF.to.longDT <- function(df_wide, return_DF = TRUE) {
  assertthat::assert_that(assertthat::is.flag(return_DF))
  assertthat::assert_that(data.table::is.data.table(df_wide)||is.data.frame(df_wide))
  # wrapping any call inside this function call will suppress all warnings:
  SuppressAllWarnings <- function(expr) {
    h <- function (w) {
      if (!is.null(w$message)) invokeRestart("muffleWarning")
    }
    withCallingHandlers(expr, warning = h)
  }
  # old attributes saved:
  newattrs <- CopyAttributes(attrslist=attributes(df_wide), dataform="long") 	# save attributes from input wide format data
  Nsamp <- nrow(df_wide)
  all_ts <- get_allts(df_wide) # calculate actual time values used in the simulated data
  attnames <- attr(df_wide, "attnames")
  # node_nms <- attr(df_wide, "node_nms")
  node_nms <- colnames(df_wide)
  # browser()
  node_nms_split <- strsplit(node_nms, "_")
  # node_nms_split2 <- strsplit(node_nms, '.', fixed = TRUE)
  node_nms_vec <- sapply(node_nms_split, '[', 1)
  # sapply(node_nms_split2, '[', 1)

  node_nms_unq <- unique(node_nms_vec)
  varying <- list()
  v.names <- NULL
  bsl.names <- NULL

  dprint("all_ts"); dprint(all_ts)
  dprint("node_nms"); dprint(node_nms)

  if (exists("setDTthreads")) setDTthreads(1)

  # if there are no time-points (t) attributes, then the long vs. wide format is undefined.
  if (length(all_ts)==0) return(df_wide)
  #******************************************************************************************************
  # this will create a copy of the object, wont modify the original object in the calling environment
  if (!is.data.table(df_wide)) df_wide <- data.table(df_wide)
  # !!!!! NOTE: THIS WILL convert the data.frame in the calling environment to data.table!!!!!!
  # To cancel that run class(df_wide) <- "data.frame" afterwards
  # data.table::setDT(df_wide, giveNames=TRUE, keep.rownames=FALSE)  # convert data.frame to data.table by reference (without copying)
  #******************************************************************************************************
  DF.to.longDT_run <- function(dat_df, t_pts, lvars, bslvars=NULL) {
    t_vec <- "_"%+%(t_pts)

    for (var_nm in lvars) { # one TV var at a time approach
      value_vars <- var_nm%+%t_vec
      SuppressAllWarnings(
        DT_melt <- melt(dat_df, id.vars = "ID", measure.vars = value_vars, variable.factor = TRUE, na.rm = FALSE)
      )
      # DT_melt <- melt.data.table(dat_df, id.vars="ID", measure.vars=value_vars, variable.factor=TRUE, na.rm=FALSE)
      # DT_melt <- data.table::melt(dat_df, id.vars="ID", measure.vars=value_vars, variable.factor=TRUE, na.rm=FALSE)
      var_nm_rep <- rep.int(var_nm, nrow(DT_melt))
      t_rep <- rep(t_pts, each=nrow(dat_df))
      DT_melt[, c("LVname","t") := list(var_nm_rep,t_rep)]
      DT_l <- data.table::dcast.data.table(DT_melt, t + ID ~ LVname)
      data.table::setkeyv(DT_l, c("ID", "t"))
      # data.table::setkey(DT_l, ID, t)
      if (which(lvars%in%var_nm)==1) {
        DT_l_fin <- DT_l
      } else {
        DT_l_fin[, (var_nm) := DT_l[[var_nm]]]
      }
    }
    if (!is.null(bslvars)) {
      DT_l_bsl <- dat_df[,unique(c("ID",bslvars)), with=FALSE]
      # data.table::setkey(DT_l_bsl, ID)  # needs to be sorted
      data.table::setkeyv(DT_l_bsl, "ID")  # needs to be sorted
      DT_l_fin <- DT_l_bsl[DT_l_fin]
    }
    DT_l_fin
  }
  # for each node, determinine if its baseline or time-varying
  # if time-varying and columns are missing for some of its time-points => impute those time-points as NA (last ifelse)
  for (cur_node_name in node_nms_unq) {
    idx_node <- sapply(node_nms_split, function(t_nodename) t_nodename[1]%in%cur_node_name)
    count_ts <- sum(idx_node)
    node_t_vals <- unlist(sapply(node_nms_split, function(t_nodename) if(t_nodename[1]%in%cur_node_name) as.integer(t_nodename[2])))
    # if node name has no "_", assume its bsl and assign first time-point:
    if (all(is.na(node_t_vals))) node_t_vals <- all_ts[1]
    dprint("cur_node_name: "%+%cur_node_name);
    dprint("all_ts"); dprint(all_ts)
    # dprint("idx_node"); dprint(idx_node)
    dprint("count_ts"); dprint(count_ts)
    dprint("node_t_vals"); dprint(node_t_vals)

    if (count_ts==length(all_ts)) { # the number of time-points is equal to all time-points =>  everything is well defined
      varying <- c(varying, list(node_nms[idx_node]))
      v.names <- c(v.names, cur_node_name)
    } else if ((count_ts==1) & (all_ts[all_ts%in%node_t_vals][1]==all_ts[1])) { # only one time-point and its the first-one
      # node is baseline: everything is fine, do nothing, just save its name
      bsl.dat.name <- node_nms[which(node_nms_vec%in%cur_node_name)][1]
      bsl.names <- c(bsl.names, cur_node_name)
      if (bsl.dat.name!=cur_node_name) {
        setnames(df_wide, bsl.dat.name, cur_node_name)
      }
    } else { # in between 1 & 2 => need to impute to length(all_ts) filling up with NAs
      ts_to_NAs <- all_ts[!all_ts%in%node_t_vals]
      ts_not_NAs <- all_ts[all_ts%in%node_t_vals]
      NA_nodenames <- cur_node_name%+%"_"%+%ts_to_NAs # define these node names as missing and add them to df
      df_wide[, (NA_nodenames) := NA]
      all_vnames <- cur_node_name%+%"_"%+%all_ts
      varying <- c(varying, list(all_vnames))
      v.names <- c(v.names, cur_node_name)
    }
  }
  dprint("v.names"); dprint(v.names)
  dprint("bsl.names"); dprint(bsl.names)
  #*** the baseline variable names are also passed as last arg
  simDT_long <- DF.to.longDT_run(dat_df=df_wide, t_pts=all_ts, lvars=v.names, bslvars=bsl.names)
  if (length(attnames)>0) v.names <- v.names[!v.names%in%attnames]
  dprint("v.names after rem attrs"); dprint(v.names)
  simDT_long <- simDT_long[rowSums(is.na(simDT_long[,v.names, with=FALSE]))<length(v.names),] # 3) remove last NA rows (after censoring or EFUP)

  # attributes(simDT_long) <- c(attributes(simDT_long), newattrs)
  newattrs_names <- names(newattrs)
  for (attr_name in newattrs_names) {
    setattr(x = simDT_long, name = attr_name, value = newattrs[[attr_name]])
  }
  # attr(simDT_long, "v.names") <- v.names
  setattr(x = simDT_long, name = "v.names", value = v.names)

  if (return_DF) class(simDT_long) <- "data.frame"

  return(simDT_long)
}

# NOT IMPLEMENTED
# will convert the dataset back to wide format from long, get the necessary attributes from simdf_long
# convert_to_wide <- function(simdf_long) {
# 	# Nsamp <- nrow(df_wide)
#     # all_ts <- get_allts(simdf_long) # calculate actual time values used in the simulated data
#     attnames <- attr(simdf_long, "attnames")
#     node_nms <- attributes(simdf_long)$names[-1]
#     # node_nms_split <- strsplit(node_nms, "_")
#     # node_nms_unq <- unique(sapply(node_nms_split, '[', 1))
# 	varying <- list()
# 	v.names <- NULL
# 	simdf_wide <- reshape(data=simdf_long, idvar=, varying=varying, v.names=v.names, timevar = "t", times = all_ts, sep="_", direction="wide")
# 	simdf_wide
# }