R/event_vector.R

In bbuchsbaum/fmrireg: R package for the anlysis of fmri data

#' @keywords internal
#' @noRd
.sanitizeName <- function(name) {
  name <- gsub(":", ".", name)
  name <- gsub(" ", "", name)
  name <- gsub("\\)$", "", name)
  name <- gsub("[\\(\\)]", "_", name, perl=TRUE)
  name <- gsub(",", "_", name)
  name <- gsub("\\.$", "", name)
  name
}

#' @keywords internal
#' @noRd
is.increasing <- function(vec) {
  all(diff(vec) >= 0)
}

#' @keywords internal
#' @noRd
is.strictly.increasing <- function(vec) {
  all(diff(vec) > 0)
}

#' @keywords internal
#' @import assertthat
#' @noRd
.checkEVArgs <- function(name, vals, onsets, blockids, durations=NULL) {

  assert_that(length(onsets) == length(vals))
  
  ## no NA onsets allowed
  assert_that(all(!is.na(onsets)))
  
  sons <- split(onsets, blockids)

  for (ons in sons) {
    assert_that(is.strictly.increasing(ons))
  }
  
  if (is.null(durations) || length(durations) == 1) {
    durations <- rep(durations, length(onsets))
  }
  
  assert_that(length(durations) == length(vals))
  assert_that(length(blockids) == length(vals))
  
  list(varname=name, value=vals, onsets=onsets, durations=durations, blockids=blockids)
}



#' Create an event model term from a named list of variables.
#'
#' This function generates an event model term from a list of named variables,
#' along with their onsets, block IDs, and durations. Optionally, a subset of
#' onsets can be retained.
#'
#' @param evlist A list of named variables.
#' @param onsets A vector of onset times for the experimental events in seconds.
#' @param blockids A vector of block numbers associated with each onset.
#' @param durations A vector of event durations (default is 1).
#' @param subset A logical vector indicating the subset of onsets to retain (default is NULL).
#' 
#' @return A list containing the following components:
#'   - varname: A character string representing the variable names, concatenated with colons.
#'   - events: A list of event variables.
#'   - subset: A logical vector indicating the retained onsets.
#'   - event_table: A tibble containing event information.
#'   - onsets: A vector of onset times.
#'   - blockids: A vector of block numbers.
#'   - durations: A vector of event durations.
#' 
#' @examples 
#' x1 <- factor(rep(letters[1:3], 10))
#' x2 <- factor(rep(1:3, each=10))
#' eterm <- event_term(list(x1=x1,x2=x2), onsets=seq(1,100,length.out=30), 
#'                     blockids=rep(1,30))
#' 
#' x1 <- rnorm(30)
#' x2 <- factor(rep(1:3, each=10))
#' eterm <- event_term(list(x1=x1,x2=x2), onsets=seq(1,100,length.out=30), 
#'                     blockids=rep(1,30), subset=x1>0)
#'
#' @export
event_term <- function(evlist, onsets, blockids, durations = 1, subset=NULL) {
  
  assert_that(is.increasing(blockids), msg="'blockids' must consist of strictly increasing integers")
            
  if (is.null(subset)) { subset=rep(TRUE, length(onsets)) }
  
  if (length(durations) == 1) {
    durations <- rep(durations, length(onsets))
  }
  
  vnames <- names(evlist)
  evs <- lapply(1:length(evlist), function(i) EV(evlist[[i]], vnames[i], 
                                                 onsets=onsets, 
                                                 blockids=blockids, 
                                                 durations=durations,
                                                 subset=subset))
  
  names(evs) <- sapply(evs, function(ev) ev$varname)

  pterms <- unlist(lapply(evs, function(ev) ev$varname))
  
  len <- sum(subset)
  

  etab <- suppressMessages(tibble::as_tibble(lapply(pterms, function(termname) {
    if (is_continuous(evs[[termname]])) {
      rep(.sanitizeName(termname), len)
    } else {
      evs[[termname]]$value
    }			
  }), .name_repair="check_unique"))
  
  names(etab) <- sapply(pterms, .sanitizeName)
  varname <- paste(sapply(evs, function(x) x$varname), collapse=":")
  
  ret <- list(varname=varname, 
              events=evs, 
              subset=subset, 
              event_table=etab, 
              onsets=evs[[1]]$onsets, 
              blockids=evs[[1]]$blockids, 
              durations=evs[[1]]$durations)
  class(ret) <- c("event_term", "event_seq")
  ret
}


#' @export
event_table.event_term <- function(x) x$event_table


#' EV
#' 
#' factory function for creating 'event' types: event_factor, event_variable, event_basis, event_matrix.
#' 
#' @param vals the event values
#' @param name the name of the event variable
#' @param onsets the event onsets.
#' @param blockids the block ids associated with each event (must be non-decreasing)
#' @param durations the duration of each event.
#' @param subset a \code{logical} vector indicating the subset of events to keep
#' 
#' @examples 
#' 
#' ev_fac <- EV(factor(c("A", "B", "C")), "fac", onsets=c(1,10,20), 
#' blockids=rep(1,3))
#' 
#' ev_fac2 <- EV(factor(c("A", "B", "C")), "fac", onsets=c(1,10,20), 
#' blockids=rep(1,3), subset=c(TRUE, TRUE, FALSE))
#' 
#' ev_numeric <- EV(c(1,2,3), "fac", onsets=c(1,10,20), 
#' blockids=rep(1,3))
#' @keywords internal
#' @noRd
EV <- function(vals, name, onsets, blockids, durations = 1, subset=rep(TRUE,length(onsets))) {
  
  if (length(durations) == 1) {
    durations <- rep(durations, length(onsets))
  }
  
  if (is.matrix(vals) && NCOL(vals) == 1) {
    vals <- vals[, 1, drop=TRUE]
  }
  
  ## subset is eagerly applied here.
  
  if (inherits(vals, "ParametricBasis")) {
    event_basis(vals, onsets, blockids, durations,subset)	
  } else if (is.factor(vals) || is.character(vals)) {
    vals <- factor(as.character(vals)[subset])
    event_factor(vals, name, onsets[subset], blockids[subset], durations[subset])
  }else if (is.numeric(vals)) {
    event_variable(vals[subset], name, onsets[subset], blockids[subset], durations[subset])
  } else if (is.matrix(vals)) {
    event_matrix(vals[subset,], name, onsets[subset], blockids[subset], durations[subset])
  } else {
    stop(paste("cannot create event_seq from type: ", typeof(vals)))
  }
  
}

#' Create a categorical event sequence from a factor
#'
#' This function generates a categorical event sequence object from a given factor.
#' It can be used to create event sequences for categorical predictors in fMRI data analyses.
#'
#' @param fac A factor representing the categorical event sequence
#' @param name A character string providing a name for the event sequence
#' @param onsets A numeric vector of onsets for each event in the sequence
#' @param blockids A numeric vector of block identifiers for each event in the sequence (default: rep(1, length(fac)))
#' @param durations A numeric vector of durations for each event in the sequence (default: rep(0, length(fac)))
#'
#' @return An object representing the categorical event sequence, with class "event_factor" and "event_seq"
#'
#' @examples 
#' efac <- event_factor(factor(c("a", "b", "c", "a", "b", "c")), "abc", onsets=seq(1, 100, length.out=6))
#'
#' @seealso \code{\link{EV}}, \code{\link{event_model}}
#' @export 
event_factor <- function(fac, name, onsets, blockids=rep(1,length(fac)), durations=rep(0, length(fac))) {
  if (!is.factor(fac)) {
    warning("argument 'fac' is not a factor, converting to factor")
    fac <- factor(as.character(fac))
  }
  
  ret <- .checkEVArgs(name, fac, onsets, blockids, durations)
  ret$continuous = FALSE
  #ret$split_onsets <- split(efac$onsets, efac$value)
  class(ret) <- c("event_factor", "event_seq")
  ret
}        

#' Create a continuous valued event sequence from a numeric vector
#'
#' This function generates a continuous valued event sequence object from a given numeric vector.
#' It can be used to create event sequences for continuous predictors in fMRI data analyses.
#'
#' @param vec A numeric vector representing the continuous event sequence values
#' @param name A character string providing a name for the event sequence
#' @param onsets A numeric vector of onsets for each event in the sequence
#' @param blockids A numeric vector of block identifiers for each event in the sequence (default: 1)
#' @param durations A numeric vector of durations for each event in the sequence (default: NULL)
#'
#' @return An object representing the continuous valued event sequence, with class "event_variable" and "event_seq"
#'
#' @examples 
#' evar <- event_variable(c(1, 2, 3, 4, 5, 6), "example_var", onsets=seq(1, 100, length.out=6))
#'
#' @seealso \code{\link{EV}}, \code{\link{event_model}}
#' @export
event_variable <- function(vec, name, onsets, blockids=1, durations=NULL) {
  stopifnot(is.vector(vec))
  
  if (is.factor(vec)) {
    stop("cannot create an event_variable from a factor, use 'event_factor'.")
  }
  
  ret <- .checkEVArgs(name, vec, onsets, blockids, durations)
  ret$continuous <- TRUE
  class(ret) <- c("event_variable", "event_seq")
  ret
  
}       

#' Create a continuous valued event set from a matrix
#'
#' This function generates a continuous valued event set object from a given matrix.
#' It is useful for creating event sequences with multiple continuous predictors for fMRI data analyses.
#'
#' @param mat A matrix of continuous event sequence values, with one row per event
#' @param name A character string providing a name for the event set
#' @param onsets A numeric vector of onsets for each event in the sequence
#' @param blockids A numeric vector of block identifiers for each event in the sequence (default: 1)
#' @param durations A numeric vector of durations for each event in the sequence (default: NULL)
#'
#' @return An object representing the continuous valued event set, with class "event_matrix" and "event_seq"
#'
#' @examples 
#' mat <- matrix(rnorm(200), 100, 2)
#' onsets <- seq(1, 1000, length.out=100)
#' durations <- rep(1, 100)
#' blockids <- rep(1, 100)
#'
#' eset <- event_matrix(mat, "eset", onsets, durations, blockids)
#'
#' @export
event_matrix <- function(mat, name, onsets, blockids=rep(1, ncol(mat)), durations=NULL) {
#  browser()
  stopifnot(is.matrix(mat))
  
  ret <- .checkEVArgs(name, as.vector(mat[,1]), onsets, blockids, durations)
  ret$continuous <- TRUE
  
  if (is.null(colnames(mat))) {
    colnames(mat) <- 1:NCOL(mat)
  }
  
  ret$value <- mat
  
  class(ret) <- c("event_matrix", "event_seq")
  ret
}


#' Create an event set from a ParametricBasis object
#'
#' This function generates an event set object from a given basis object of type `ParametricBasis`.
#' It is useful for creating event sequences based on basis functions for fMRI data analyses.
#'
#' @param basis A ParametricBasis object containing the basis functions
#' @param onsets A numeric vector of onsets for each event in the sequence
#' @param blockids A numeric vector of block identifiers for each event in the sequence (default: 1)
#' @param durations A numeric vector of durations for each event in the sequence (default: NULL)
#' @param subset A logical vector indicating a subset of the basis object to use (default: TRUE for all elements)
#'
#' @return An object representing the event set based on the basis functions, with class "event_basis" and "event_seq"
#'
#' @import assertthat
#' @examples 
#' # Create a ParametricBasis object
#' basis <- ParametricBasis("Gamma", shape = 6, rate = 0.9)
#' onsets <- seq(0, 20, length.out = 21)
#' blockids <- rep(1, length(onsets))
#'
#' # Generate an event_basis object
#' ebasis <- event_basis(basis, onsets, blockids)
#'
#' @export
event_basis <- function(basis, onsets, blockids=1, durations=NULL, subset=rep(TRUE, length(onsets))) {
  assertthat::assert_that(inherits(basis, "ParametricBasis"))
  
  if (any(!subset)) {
    basis <- sub_basis(basis, subset)
  }

  ret <- .checkEVArgs(basis$name, basis$y[,1], onsets[subset], blockids[subset], durations[subset])
  #basis$y <- basis$y[subset,]
  ret$value <- basis$y
  ret$continuous <- TRUE
  ret$basis <- basis
  class(ret) <- c("event_basis", "event_seq")
  ret
}


#' @export
#' @rdname levels
levels.event_factor <- function(x) levels(x$value) 

#' @export
#' @rdname levels
levels.event_variable <- function(x) x$varname 

#' @export
#' @rdname levels
levels.event_matrix <- function(x) colnames(x$value) 

#' @export
#' @rdname levels
levels.event_set <- function(x) colnames(x$value) 

#' @export
levels.event_basis <- function(x) levels(x$basis)

#' @export
formula.event_term <- function(x, ...) as.formula(paste("~ ", "(", paste(parent_terms(x), collapse=":"), 
                                                   "-1", ")"))

#' @export
#' @rdname levels
levels.event_term <- function(x) {
  facs <- x$events[!sapply(x$events, is_continuous)]
  if (length(facs) == 1) {
    levels(facs[[1]])
  } else {
    facs <- lapply(facs, function(x) x$value)
    f <- function(...) interaction(..., drop=TRUE, sep=":")
    levels(do.call(f, facs))
  }
}


#' Retrieve cells of an event_factor object.
#'
#' This function extracts cells from an event_factor object, optionally removing
#' empty cells. Note that the removal of empty cells is not implemented.
#'
#' @param x An event_factor object.
#' @param drop.empty Logical. If TRUE, empty cells will be removed (not implemented).
#' @param ... Additional arguments to be passed to the function.
#'
#' @return A list of data frames containing the cells of the event_factor object.
#'
#' @export
#' @rdname cells
cells.event_factor <- function(x, drop.empty=TRUE,...) {
  etab <- data.frame(onsets=x$onsets, durations=x$durations, blockids=x$blockids)
  split(etab, x$value)
}


#' Retrieve cells of an event_term object.
#'
#' This function extracts cells from an event_term object, optionally removing
#' empty cells.
#'
#' @param x An event_term object.
#' @param drop.empty Logical. If TRUE, empty cells will be removed.
#' @param ... Additional arguments to be passed to the function.
#'
#' @return A tibble containing the cells of the event_term object, with an
#'   additional "count" attribute.
#'
#' @export
#' @rdname cells
#' @examples
#' 
#' evlist <- list(fac1=factor(c("A", "B", "A", "B")), 
#'                fac2=factor(c("1", "1", "2", "2")))
#' eterm <- event_term(evlist,onsets=1:4, blockids=rep(1,4))
#' cells(eterm)
cells.event_term <- function(x, drop.empty=TRUE,...) {
  evtab <- x$event_table
  evset <- suppressMessages(tibble::as_tibble(expand.grid(lapply(x$events, levels)), .name_repair="check_unique"))
  
  which.cat <- which(!sapply(x$events, is_continuous))
  
  if (length(which.cat) > 0) {
    evs <- suppressMessages(tibble::as_tibble(lapply(evset[,which.cat], as.character), .name_repair="check_unique"))
    evt <- suppressMessages(tibble::as_tibble(lapply(evtab[,which.cat], as.character), .name_repair="check_unique"))
  
    counts <- apply(evs, 1, function(row1) {
      sum(apply(evt, 1, function(row2) {										
        all(as.character(row1) == as.character(row2))
      }))
    })
  
    if (drop.empty) {
      evset <- evset[counts > 0,,drop=F]
      attr(evset, "count") <- counts[counts > 0]
    } else {
      attr(evset, "count") <- counts
    }
  } else {
    attr(evset, "count") <- nrow(evtab)
  }
  
  evset
}


#' @noRd
#' @keywords internal
.event_set <- function(x, exclude_basis=FALSE) {
  evtab <- event_table(x)
  
  evset <- if (nbasis(x) > 1 & !exclude_basis) {
    evlist <- c(list(factor(paste("basis", 1:nbasis(x), sep = ""))), cells(x$evterm))
    names(evlist) <- c("basis", parent_terms(x$evterm))
    evlist <- lapply(evlist, levels)
    ret <- expand.grid(evlist, stringsAsFactors = TRUE)
    ret[c(2:length(ret), 1)]
  } else {
    cells(x$evterm)
  }
  
}

#' @export
#' @rdname cells
cells.covariate_convolved_term <- function(x,...) {
  unique(event_table(x))
}

#' @export
#' @importFrom stringr str_trim
cells.convolved_term <- function(x, exclude_basis=FALSE,...) {
  evtab <- event_table(x)
  evset <- .event_set(x, exclude_basis=exclude_basis)
  
  strels <- apply(apply(evtab, 2, stringr::str_trim), 1, paste, collapse = ":")
  
  strlevs <- if (nrow(evset) > 1) {
    apply(apply(evset, 2, stringr::str_trim), 1, paste, collapse = ":")
  } else {
    as.character(evset[1,1])
  }
  
  attr(evset, "rownames") <- strlevs
  
  counts <- if (exclude_basis) {
    rep(attr(cells(x$evterm), "count"), each = 1)
  } else {
    rep(attr(cells(x$evterm), "count"), each = nbasis(x))
  }
  
  ret <- evset[counts > 0, , drop = F]
  attr(ret, "count") <- counts[counts > 0]
  ret
  
}

#' @export
#' @rdname conditions
conditions.fmri_term <- function(x, ...) {
  colnames(design_matrix(x))
}

#' @export
#' @rdname conditions
conditions.convolved_term <- function(x,...) {
  colnames(design_matrix(x))
}

#' @export
#' @rdname conditions
conditions.afni_hrf_convolved_term <- function(x,...) {
  conditions(x$evterm)
}

#' @export
#' @rdname conditions
conditions.afni_trialwise_convolved_term <- function(x,...) {
  conditions(x$evterm)
}


#' @export
#' @rdname conditions
conditions.event_term <- function(x, drop.empty=TRUE, ...) {
  
  .cells <- cells(x, drop.empty=drop.empty)
  pterms <- parent_terms(x)
  levs <- apply(.cells, 1, paste, collapse=":")
  
  splitlevs <- strsplit(levs, ":")
  ret <- lapply(1:length(pterms), function(i) {
    lev <- sapply(splitlevs, "[[", i)
    term <- pterms[[i]]
    
    if (length(levels(x$events[[i]])) > 1) {
      paste(.sanitizeName(pterms[i]), "[", lev, "]", sep="")
    } else {
      .sanitizeName(pterms[i])
    }
  })
  
  do.call(function(...) paste(..., sep=":"), ret)
}


## TODO is columns ever used? do we need this function

#' @export
#' @rdname columns
columns.event_term <- function(x) as.vector(unlist(lapply(x$events, columns)))

#' @export
#' @rdname columns
columns.event_seq <- function(x) x$varname

#' @export
#' @rdname columns
columns.event_matrix <- function(x) paste0(.sanitizeName(x$varname), ".", levels(x))

#' @export
#' @rdname columns
columns.event_set <- function(x) paste0(.sanitizeName(x$varname), ".", levels(x))

#' @export
#' @rdname columns
columns.event_basis <- function(x) columns(x$basis)


#' @export
#' @rdname parent_terms
parent_terms.event_term <- function(x) unlist(lapply(x$events, function(ev) ev$varname))


#' @export
is_continuous.event_seq <- function(x) x$continuous

#' @export
is_continuous.event_factor <- function(x) x$continuous

#' @export
is_categorical.event_seq <- function(x) !x$continuous

#' @export
is_continuous.event_term <- function(x) all(sapply(x$events, function(x) is_continuous(x)))

#' @export
is_categorical.event_term <- function(x) !is_continuous(x)


#' @export
is_categorical.event_seq <- function(x) !x$continuous


#' @export
elements.event_matrix <- function(x, values=TRUE, ...) {
  if (values) {
    ret <- x$value
    colnames(ret) <- colnames(x)
    ret <- list(ret)
    names(ret) <- .sanitizeName(x$varname)
    ret
  } else {
    N <- length(x$onsets)
    vnames <- colnames(x)
    res <- lapply(vnames, function(el) rep(el, N))
    mat <- do.call(cbind, res)
    colnames(mat) <- vnames			
    ret <- list(mat)
    names(ret) <- .sanitizeName(x$varname)
    ret			
  }
}


#' @export
elements.event_seq <- function(x, values = TRUE, ...) {
  if (values) {
    ret <- list(x$value)
    names(ret) <- x$varname
    ret
  } else {
    #ret <- list(rep(varname(x), length(x)))
    ret <- list(rep(x$varname, length(x)))
    names(ret) <- x$varname
    ret
  }
}

#' @export
elements.event_basis <- function(x, values=TRUE, transformed=TRUE, ...) {
  if (values && !transformed) {
    x$value$x				
  } else if (values) {
    #browser()
    ret <- x$basis$y
    colnames(ret) <- columns(x)
    n <- .sanitizeName(x$varname)
    ret <- list(ret)
    names(ret) <- n
    ret
  } else {
    N <- length(x)
    vnames <- columns(x)
    res <- lapply(vnames, function(el) rep(el, N))
    mat <- do.call(cbind, res)
    colnames(mat) <- vnames			
    ret <- list(mat)
    names(ret) <- .sanitizeName(x$varname)
    ret		
  }
}


#' @export
elements.event_term <- function(x, values=TRUE, ...) {
  els <- lapply(x$events, elements, values=values)
  n <- sapply(names(els), function(nam) .sanitizeName(nam))
  names(els) <- as.vector(n)
  els
}

#' @export
onsets.convolved_term <- function(x) {
  onsets(x$evterm)
}

#' @export
onsets.event_term <- function(x) {
  x$onsets
}

#' @export
blockids.event_term <- function(x) {
  x$blockids
}

#' @export
blockids.convolved_term <- function(x) {
  blockids(x$evterm)
}

#' Split onsets of an event_term object
#'
#' This function splits the onsets of an event_term object based on its factor levels or block identifiers.
#' It is useful for processing fMRI data when analyzing event-related designs with multiple conditions or blocks.
#'
#' @param x An event_term object
#' @param sframe A data frame representing the sampling frame
#' @param global A logical value indicating whether to use global onsets (default: FALSE)
#' @param blocksplit A logical value indicating whether to split onsets by block identifiers (default: FALSE)
#' @param ... Additional arguments passed to other methods
#'
#' @return A list of numeric vectors representing the split onsets for each factor level or block identifier
#'
#' @export
split_onsets.event_term <- function(x, sframe, global=FALSE,blocksplit=FALSE, ...) {
  ### need to check for 0 factors
  facs <- x$events[!sapply(x$events, is_continuous)]
  
  if (length(facs) == 0) {
    ons <- if (global) {
      global_onsets(sframe, onsets(x), blockids(x))
    } else {
      onsets(x)
    }
    return(list(split(ons, blockids(x))))
    
  }
  
  facs <- lapply(facs, function(fac) unlist(elements(fac)))
  
  f <- function(...) {
    interaction(..., drop=TRUE, sep=":")
  }
            
  cfac <- try(do.call(f, facs))
  #if (inherits(cfac, "try-error")) {
  #  stop("could not construct crossed factor for event_term `x`")
  #}
            
  ret <- if (global) {
    split(global_onsets(sframe, onsets(x), blockids(x)), cfac)
  } else {
    split(onsets(x), cfac)
  }
  
  if (blocksplit) {
    bsplit <- split(blockids(x), cfac)
    ret <- lapply(1:length(ret), function(i) {
      split(ret[[i]], bsplit[[i]])
    })
  }
  
  names(ret) <- longnames(x)
  ret
}




#' Convolve HRF with Design Matrix
#'
#' This function convolves a given hemodynamic response function (HRF) with the design matrix of an fMRI study.
#' It is useful for modeling the expected BOLD signal in response to experimental conditions.
#'
#' @param hrf A numeric vector representing the hemodynamic response function
#' @param dmat A design matrix with columns representing different experimental conditions
#' @param globons A numeric vector of global onsets for each event
#' @param durations A numeric vector of event durations
#' @param summate A logical value indicating whether to summate the convolved HRF (default: TRUE)
#'
#' @return A list of regressors, one for each column in the design matrix.
#' @export
convolve_design <- function(hrf, dmat, globons, durations, summate=TRUE) {
  cond.names <- names(dmat)
  
  if (any(is.na(dmat)) || any(is.na(globons))) {
    keep <- apply(dmat, 1, function(vals) all(!is.na(vals)))
    keep[is.na(globons)] <- FALSE
    dmat <- dmat[keep,]
    durations <- durations[keep]
    globons <- globons[keep]
  } 
  
  reglist <- purrr::map(1:ncol(dmat), function(i) {
    amp <- dmat[,i][[1]]
    nonzero <- which(amp != 0)
    ## issue with single_trial_regressor
    if (length(nonzero) == 0) {
      null_regressor(hrf)
      ## scaling issue with single_trial_regressor
      #} #else if (length(nonzero) == 1) {
      #single_trial_regressor(globons[nonzero], hrf, amplitude=amp[nonzero], duration=durations[nonzero])
    } else {
      regressor(globons[nonzero], hrf, amplitude=amp[nonzero], duration=durations[nonzero], summate=summate)
    }
  })
  
  reglist
  
}

#' Convolve an event-related design matrix with an HRF.
#'
#' This function takes an event-related design matrix and convolves it with
#' a specified Hemodynamic Response Function (HRF) to create a new design matrix
#' suitable for fMRI analysis. It also supports additional arguments for
#' flexibility and customization.
#'
#' @importFrom tibble as_tibble
#' @importFrom magrittr %>%
#' @importFrom dplyr group_by select do ungroup
#' @autoglobal
#' @export
#'
#' @param x A data frame containing the input design matrix.
#' @param hrf A Hemodynamic Response Function to convolve the design matrix with.
#' @param sampling_frame A data frame specifying the sampling frame for the analysis.
#' @param drop.empty Logical. If TRUE, empty rows in the design matrix will be removed.
#' @param summate Logical. If TRUE, the convolved design matrix will be summed.
#' @param precision Numeric. The desired precision for the calculations.
#' @param ... Additional arguments to be passed to the function.
#'
#' @return A convolved design matrix, in tibble format.
convolve.event_term <- function(x, hrf, sampling_frame, drop.empty=TRUE, 
                                summate=TRUE, precision=.3,...) {
  
  
  globons <- global_onsets(sampling_frame, x$onsets, x$blockids)
  
  durations <- x$durations
  blockids <- x$blockids
  
  nimages <- sum(sampling_frame$blocklens)
  
  cnames <- conditions(x)
  dmat <- design_matrix(x, drop.empty)
  ncond <- ncol(dmat)

  #print("convolving")
 
  cmat <- dmat |> dplyr::mutate(.blockids=blockids, .globons=globons, .durations=durations) |> 
    dplyr::group_by(.blockids) %>%
    dplyr::do({
      d <- dplyr::select(., 1:ncond)
      reg <- convolve_design(hrf, d, .$.globons, .$.durations, summate=summate)
      sam <- samples(sampling_frame, blockids=as.integer(as.character(.$.blockids[1])), global=TRUE)
      
      ## TODO could bee parallelized
      ## ret <- do.call(rbind, furrr::future_map(reg, function(r) evaluate(r, sam))) 
      ret <- do.call(cbind, lapply(seq_along(reg), function(ri) {
        vname <- paste0("v", ri)
        evaluate(reg[[ri]], sam, precision=precision) 
        
        #%>% select({{vname}} := value)
        #names(tmp) <- paste0("v", ri)
        #suppressMessages(as_tibble(tmp))
      })) 
      
      ret <- suppressMessages(tibble::as_tibble(ret, .name_repair="minimal"))
      names(ret) <- paste0("v", 1:length(reg))
      ret
  }) 
  
  cmat <- cmat %>% dplyr::ungroup() %>% dplyr::select(-.blockids)
  
  if (nbasis(hrf) > 1) {
    blevs <- paste("[", 1:nbasis(hrf), "]", sep="")
    cnames <- unlist(lapply(cnames, function(prefix) paste(prefix, ":basis", blevs, sep="")))
  } 
            
  
  colnames(cmat) <- cnames
  suppressMessages(tibble::as_tibble(cmat, .name_repair="check_unique"))
  
  
}

#' Compute F-contrasts for Event Term
#'
#' This function computes F-contrasts for an event term object, considering main effects and interactions.
#'
#' @param x An event term object
#' @param ... Additional arguments passed to the function
#'
#' @return A list of contrast matrices for main effects and interactions
#' @export
Fcontrasts.event_term <- function(x,...) {
  cellcount <- attr(cells(x, drop.empty=FALSE), "count")
  if (any(cellcount) == 0) {
    stop("currently cannot compute Fcontrasts for non-orthogonal design.")
  }
  ##browser()
  ## TODO check for no empty cells, otherwise everything fails
  which_cat <- which(sapply(x$events, function(obj) is_categorical(obj)))
  assert_that(length(which_cat) > 0, msg="Fcontrasts cannot be computed for terms with no categorical variables")
  ## factors comprising this term
  pterms <- parent_terms(x)[which_cat]
  evs <- x$events[which_cat]
  cond <- conditions(x)
  facnames <- names(evs)
  
  Clist <- lapply(evs, function(ev) rep(1, length(levels(ev))))
  Dlist <- lapply(evs, function(ev) t(-diff(diag(length(levels(ev))))))
  
  nfac <- length(Clist)
  
  valid_cells <- cellcount > 0

  main_effects <- lapply(length(Clist):1, function(i) {
    #print(i)
    Dcon <- Dlist[[i]]
    Cs <- Clist[-i]
    mats <- vector(nfac, mode="list")
    mats[[i]] <- Dcon
    mats[seq(1, nfac)[-i]] <- Cs
    ret <- Reduce(kronecker, rev(mats))
    if (!(all(valid_cells))) {
      ret <- ret[valid_cells,,drop=FALSE]
      if (ncol(ret) > 1) {
        ret <- svd(ret)$u
      } else {
        ret <- scale(ret, center=TRUE, scale=FALSE)
      }
    }
    row.names(ret) <- cond
    ret
  })
  
  names(main_effects) <- rev(facnames)
  
  if (length(facnames) > 1 && all(valid_cells)) {
    interactions <- vector(length(Clist)-1, mode="list")
    for (i in length(Clist):2) {
      icomb <- combn(nfac, i)
      ret <- lapply(1:ncol(icomb), function(j) {
        ind <- icomb[,j]
        mats <- vector(nfac, mode="list")
        mats[ind] <- Dlist[ind]
        if (length(ind) < nfac) {
          mats[-ind] <- Clist[-ind]
        }
      
        cmat <- Reduce(kronecker, mats)
        row.names(cmat) <- cond
        cmat
      })
      cnames <- apply(icomb, 2, function(i) paste0(facnames[i], collapse=":"))
      names(ret) <- cnames
      interactions[[i-1]] <- ret
    }
    
    return(c(main_effects, unlist(interactions, recursive=FALSE)))
  } else {
    main_effects
  }
  
}
                 

#' @importFrom tibble as_tibble
#' @importFrom purrr map_chr
#' @export
design_matrix.event_term <- function(x, drop.empty=TRUE,...) {

  locenv <- new.env()
  pterms <- map_chr(parent_terms(x), .sanitizeName)	
  
  for (ev in x$events) {
    vname <- .sanitizeName(ev$varname)
    els <- elements(ev, values=TRUE)
    #if (length(els) == 1 && is.matrix(els[[1]])) {
    ## TODO remove special casing for matrix elements 
    #  mat <- els[[1]]
    #  for (i in 1:ncol(mat)) {
    #    assign(columns(ev)[i], mat[,i],envir=locenv)
    #  }
    #} else {
      lapply(names(els), function(n) assign(n, els[[n]],envir=locenv))
    #}
  }
  
  
  els <- as.data.frame(elements(x))
  nas <- try(apply(els,1, function(vals) any(is.na(vals))))
  counts <- attr(cells(x, drop=FALSE), "count")
  
  #print(ncol(els))
  #browser()
  mat <- if (ncol(els) == 1 && is.factor(els[,1]) && length(levels(els[,1])) == 1) {
    ## a 1 level term
    cbind(rep(1, NROW(els))) 
  } else { 		
    model.matrix(formula(x), data=locenv) 
  }
  

  ### multiply design matrix by subset
  rmat <- mat #* x$subset
  
  #remove rows with NAS
  if (any(nas)) {
    rmat <- matrix(0, nrow(x$event_table), length(conditions(x, drop.empty=FALSE)))
    rmat[!nas,] <- mat
    rmat[nas,] <- NA				
  } 
  

  # remove columns with no events (postpone this to later stage?) 
  if (any(counts == 0) && (length(conditions(x, drop=F)) == length(counts)) && drop.empty) {
    rmat <- rmat[, !(counts==0), drop=FALSE]
    colnames(rmat) <- conditions(x, drop=T)
  } else {
    colnames(rmat) <- conditions(x, drop=F)			
  }
  
  suppressMessages(tibble::as_tibble(rmat, .name_repair="check_unique"))
}



#' @export
print.event_term <- function(x, ...) {
  cat("event_term", "\n")
  cat("  ", "Term Name: ", x$varname, "\n")
  cat("  ", "Formula:  ", as.character(formula(x)), "\n")
  cat("  ", "Num Events: ", nrow(x$event_table), "\n")
  cat("  ", "Term Types: ", paste(map_chr(x$events, ~ class(.)[[1]])))
  cat("\n")
}

#' @export
print.fmri_term <- function(x,...) {
  cat("fmri_term: ", class(x)[[1]], "\n")
  cat("  ", "Term Name: ", x$varname, "\n")
  cat("  ", "Num Rows: ", nrow(design_matrix(x)), "\n")
  cat("  ", "Num Columns: ", ncol(design_matrix(x)), "\n")
}

#' @export
print.convolved_term <- function(x,...) {
  cat("fmri_term: ", class(x)[[1]], "\n")
  cat("  ", "Term Name: ", x$varname, "\n")
  cat("  ", "Formula:  ", as.character(formula(x$evterm)), "\n")
  cat("  ", "Num Events: ", nrow(x$evterm$event_table), "\n")
  cat("  ", "Num Rows: ", nrow(design_matrix(x)), "\n")
  cat("  ", "Num Columns: ", ncol(design_matrix(x)), "\n")
  cat("  ", "Conditions: ", conditions(x), "\n")
  cat("  ", "Term Types: ", paste(map_chr(x$evterm$events, ~ class(.)[[1]])))
  cat("\n")
}

#' @export
print.afni_hrf_convolved_term <- function(x,...) {
  cat("fmri_term: ", class(x)[[1]], "\n")
  cat("  ", "Term Name: ", x$varname, "\n")
  cat("  ", "Formula:  ", as.character(formula(x$evterm)), "\n")
  cat("  ", "Num Events: ", nrow(x$evterm$event_table), "\n")
  cat("  ", "Conditions: ", conditions(x), "\n")
  cat("  ", "Term Types: ", paste(map_chr(x$evterm$events, ~ class(.)[[1]])))
  cat("\n")
}