diff --git a/.Rbuildignore b/.Rbuildignore index fd65249..c541066 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -15,3 +15,4 @@ ^Dockerfile$ ^Apptainer\.def$ ^README\.Rmd$ +^README\.html$ diff --git a/.circleci/config.yml b/.circleci/config.yml index 7cf8aac..e28f033 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -11,10 +11,24 @@ commands: r-cran-devtools \ r-bioc-rhdf5 \ r-bioc-delayedarray \ + r-cran-tiledb \ + r-cran-jsonlite \ pandoc - run: name: Install package dependencies command: R -e "devtools::install_deps(dep = TRUE, dependencies = TRUE)" + - run: + # TileDBArray is Bioconductor-only, so install_deps() cannot resolve + # it from CRAN. Without this the TileDB tests skip silently and the + # backend reads as covered when it is not. + name: Install TileDBArray (Bioconductor) + command: | + R -e 'if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")' + R -e 'BiocManager::install("TileDBArray", ask = FALSE, update = FALSE)' + - run: + name: Verify optional TileDB stack is installed + command: | + R -q -e 'for (p in c("tiledb", "TileDBArray", "jsonlite")) if (!requireNamespace(p, quietly = TRUE)) stop("TileDB test dependency missing: ", p)' jobs: check_package: diff --git a/.gitignore b/.gitignore index 64186f6..1c7bd54 100644 --- a/.gitignore +++ b/.gitignore @@ -143,6 +143,9 @@ ReadingMaterial/ vignettes/*.html vignettes/*.pdf +# knitted README (README.md is the committed artifact) +/README.html + # OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3 .httr-oauth diff --git a/Apptainer.def b/Apptainer.def index 44a6d06..5dc79ae 100644 --- a/Apptainer.def +++ b/Apptainer.def @@ -1,5 +1,5 @@ Bootstrap: docker -From: rocker/r2u:jammy +From: rocker/r2u:noble %help ModelArray - an R package for statistical analysis of fixel-wise data and beyond @@ -12,10 +12,11 @@ From: rocker/r2u:jammy org.label-schema.schema-version "1.0" %environment - export DEBIAN_FRONTEND=noninteractive + export PATH="/opt/modelarrayio/bin:${PATH}" %post set -e + sed -i 's#https://cloud.r-project.org#https://cran.r-project.org#' /etc/apt/sources.list.d/cran.sources apt-get update \ && apt-get install -y --no-install-recommends \ r-bioc-delayedarray \ @@ -37,12 +38,20 @@ From: rocker/r2u:jammy r-cran-tibble \ r-cran-tidyr \ r-cran-tidyverse \ + git \ + python3-venv \ && apt-get clean \ - && echo 'options(bspm.sudo = TRUE)' >> /etc/R/Rprofile.site \ && rm -rf /var/lib/apt/lists/* + R -e 'install.packages(c("jsonlite", "tiledb", "BiocManager"))' + R -e 'BiocManager::install("TileDBArray", ask = FALSE, update = FALSE)' + python3 -m venv /opt/modelarrayio + /opt/modelarrayio/bin/pip install --no-cache-dir git+https://github.com/PennLINC/ModelArrayIO.git + cd /ModelArray R -e 'devtools::install()' + R -e 'library(ModelArray); stopifnot(requireNamespace("tiledb", quietly = TRUE)); stopifnot(requireNamespace("TileDBArray", quietly = TRUE))' + /opt/modelarrayio/bin/modelarrayio --version %files . /ModelArray @@ -55,5 +64,3 @@ From: rocker/r2u:jammy else exec R fi - - diff --git a/DESCRIPTION b/DESCRIPTION index b6bea25..7c414aa 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -42,12 +42,15 @@ Imports: RoxygenNote: 7.3.3 Roxygen: list(markdown = TRUE) Suggests: + jsonlite, rmarkdown, knitr, testthat (>= 3.0.0), styler (>= 1.6.2), lintr, - stringr + stringr, + tiledb, + TileDBArray Config/testthat/edition: 3 VignetteBuilder: knitr URL: https://pennlinc.github.io/ModelArray, https://pennlinc.github.io/ModelArray/ diff --git a/Dockerfile b/Dockerfile index 3bd9aee..4be677e 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,11 +1,42 @@ -FROM rocker/r2u:jammy +FROM rocker/r2u:noble # Install tricky bioconductor packages and minimal LaTeX for PDF generation -RUN apt update \ - && apt install -y --no-install-recommends \ - r-cran-devtools \ - r-bioc-rhdf5 \ - r-bioc-delayedarray +RUN sed -i 's#https://cloud.r-project.org#https://cran.r-project.org#' /etc/apt/sources.list.d/cran.sources \ + && apt-get update \ + && apt-get install -y --no-install-recommends \ + r-bioc-delayedarray \ + r-bioc-hdf5array \ + r-cran-broom \ + r-cran-crayon \ + r-cran-devtools \ + r-cran-doparallel \ + r-cran-dplyr \ + r-cran-glue \ + r-cran-gratia \ + r-cran-hdf5r \ + r-cran-hdf5r.extra \ + r-cran-magrittr \ + r-cran-mgcv \ + r-cran-pbapply \ + r-cran-pbmcapply \ + r-bioc-rhdf5 \ + r-cran-tibble \ + r-cran-tidyr \ + r-cran-tidyverse \ + git \ + python3-venv \ + && apt-get clean \ + && rm -rf /var/lib/apt/lists/* + +## Install optional TileDB backend dependencies for ModelArray +RUN R -e 'install.packages(c("jsonlite", "tiledb", "BiocManager"))' \ + && R -e 'BiocManager::install("TileDBArray", ask = FALSE, update = FALSE)' + +## Install ModelArrayIO (Python package and modelarrayio CLI) +RUN python3 -m venv /opt/modelarrayio \ + && /opt/modelarrayio/bin/pip install --no-cache-dir git+https://github.com/PennLINC/ModelArrayIO.git + +ENV PATH="/opt/modelarrayio/bin:${PATH}" ## Install ModelArray (R package) @@ -13,6 +44,10 @@ COPY . /ModelArray WORKDIR /ModelArray RUN R -e 'devtools::install()' +## Verify the R and Python TileDB paths are available +RUN R -e 'library(ModelArray); stopifnot(requireNamespace("tiledb", quietly = TRUE)); stopifnot(requireNamespace("TileDBArray", quietly = TRUE))' \ + && /opt/modelarrayio/bin/modelarrayio --version + ## Add metadata: ARG BUILD_DATE ARG VCS_REF @@ -30,4 +65,4 @@ LABEL org.label-schema.build-date=$BUILD_DATE \ # but someone says it is "git branch name"?? ref: https://guide.opencord.org/cord-5.0/build_images.html org.label-schema.schema-version="1.0" # ^^these information can be viewed by: - # docker inspect pennlinc/modelarray_confixel: \ No newline at end of file + # docker inspect pennlinc/modelarray_confixel: diff --git a/NAMESPACE b/NAMESPACE index d72a0c6..28857fa 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,8 +1,10 @@ # Generated by roxygen2: do not edit by hand +S3method(print,ModelArraySummary) S3method(print,h5summary) export("%>%") export(ModelArray) +export(ModelArraySummary) export(ModelArray.gam) export(ModelArray.lm) export(ModelArray.wrap) diff --git a/R/ModelArray-package.R b/R/ModelArray-package.R index 508ac21..fc39b24 100644 --- a/R/ModelArray-package.R +++ b/R/ModelArray-package.R @@ -3,12 +3,13 @@ #' @description #' The ModelArray package provides an S4 class and associated methods for #' performing massively univariate statistical analyses on element-wise -#' (fixel, voxel, or vertex) neuroimaging data stored in HDF5 files. +#' (fixel, voxel, or vertex) neuroimaging data stored in HDF5 files or +#' TileDB stores. #' #' @details #' The core workflow is: #' \enumerate{ -#' \item Inspect an HDF5 file with \code{\link{h5summary}()} +#' \item Inspect storage with \code{\link{ModelArraySummary}()} #' \item Load data with \code{\link{ModelArray}()} #' \item Fit models with \code{\link{ModelArray.lm}}, #' \code{\link{ModelArray.gam}}, or \code{\link{ModelArray.wrap}} @@ -28,7 +29,7 @@ #' \linkS4class{ModelArray}, \code{\link{ModelArray}}, #' \code{\link{ModelArray.lm}}, \code{\link{ModelArray.gam}}, #' \code{\link{ModelArray.wrap}}, \code{\link{mergeModelArrays}}, -#' \code{\link{h5summary}} +#' \code{\link{ModelArraySummary}} #' #' @section Centralized imports: #' The following imports are consolidated here because they cannot be diff --git a/R/ModelArray_Constructor.R b/R/ModelArray_Constructor.R index 0cc69ac..b68492e 100644 --- a/R/ModelArray_Constructor.R +++ b/R/ModelArray_Constructor.R @@ -2,19 +2,19 @@ #' ModelArray class #' #' ModelArray is an S4 class that represents element-wise scalar data and -#' associated statistical results backed by an HDF5 file on disk. +#' associated statistical results backed by an HDF5 file or TileDB store on disk. #' #' @description #' A ModelArray wraps one or more element-wise scalar matrices (e.g., FD, FC, #' log_FC for fixel data) read lazily via \pkg{DelayedArray}, along with any #' previously saved analysis results. The object holds references to the -#' underlying HDF5 file and reads data on demand, making it suitable for +#' underlying storage and reads data on demand, making it suitable for #' large-scale neuroimaging datasets. #' #' @details -#' Each scalar in the HDF5 file is stored at \code{/scalars//values} +#' Each scalar is stored at \code{/scalars//values} #' as a matrix of elements (rows) by source files (columns). Source filenames -#' are read from HDF5 attributes or companion datasets. Analysis results, if +#' are read from storage metadata or companion datasets. Analysis results, if #' present, live under \code{/results//results_matrix}. #' #' ModelArray objects are typically created with the \code{\link{ModelArray}} @@ -29,7 +29,9 @@ #' @slot results A named list of analysis results. Each element is itself a #' list containing at minimum \code{results_matrix} (a #' [DelayedArray::DelayedArray][DelayedArray-class]). -#' @slot path Character. Path(s) to the HDF5 file(s) on disk. +#' @slot path Character. Path(s) to the HDF5 file(s) or TileDB store(s) on disk. +#' @slot backend Character. Resolved storage backend(s) for \code{path}: +#' \code{"hdf5"} or \code{"tiledb"}. #' #' @seealso \code{\link{ModelArray}} for the constructor, #' \code{\link{ModelArray.lm}}, \code{\link{ModelArray.gam}}, @@ -47,41 +49,49 @@ ModelArray <- setClass( results = "list", sources = "list", scalars = "list", - path = "character" - ) + path = "character", + backend = "character" + ), + prototype = list(backend = "auto") ) #' ModelArraySeed #' -#' Generates a "seed" for the h5 file format. A wrapper around HDF5ArraySeed -#' used to instantiate a delayed array +#' Generates a DelayedArray seed for the requested storage backend. #' -#' @param filepath Path to an existing h5 file. -#' @param name Name of the group/field in the h5 file. +#' @param filepath Path to an existing h5 file or TileDB store. +#' @param name Relative path to the dataset/array. #' @param type Type of DelayedArray object, used as an argument for `HDF5Array::HDF5ArraySeed`. +#' @param backend Storage backend. #' @noRd #' -ModelArraySeed <- function(filepath, name, type = NA) { +ModelArraySeed <- function(filepath, name, type = NA, backend = c("hdf5", "tiledb")) { # NOTE: the checker for if h5 groups fixels/voxels/scalars exist # (a.k.a valid fixel-wise data) is deleted, as ModelArray is generalized to any modality. - seed <- HDF5Array::HDF5ArraySeed(filepath, name = name, type = type) # HDF5Array is also from BioConductor... + backend <- match.arg(backend) + if (identical(backend, "hdf5")) { + seed <- HDF5Array::HDF5ArraySeed(filepath, name = name, type = type) # HDF5Array is also from BioConductor... + } else { + .require_tiledb_support() + seed <- TileDBArray::TileDBArraySeed(.modelarray_path(filepath, name), attr = "values") + } seed } -#' Load element-wise data from an HDF5 file +#' Load element-wise data from an HDF5 file or TileDB store #' #' @description #' Reads scalar matrices and (optionally) saved analysis results from -#' an HDF5 file and returns a \linkS4class{ModelArray} object. +#' an HDF5 file or TileDB store and returns a \linkS4class{ModelArray} object. #' #' @details #' The constructor reads each scalar listed in \code{scalar_types} from #' \code{/scalars//values}, wrapping them as #' [DelayedArray::DelayedArray][DelayedArray-class] objects. Source filenames are extracted -#' from HDF5 attributes or companion datasets. +#' from storage metadata or companion datasets. #' #' If \code{analysis_names} is non-empty, saved results are loaded from #' \code{/results//results_matrix}. @@ -92,17 +102,23 @@ ModelArraySeed <- function(filepath, name, type = NA) { #' \code{rhdf5::h5ls(filepath)}. #' #' @param filepath Character. Path to an existing HDF5 (\code{.h5}) -#' file containing element-wise scalar data. +#' file or TileDB (\code{.tdb}) store containing element-wise scalar data. #' @param scalar_types Character vector. Names of scalar groups to read -#' from \code{/scalars/} in the HDF5 file. Default is \code{c("FD")}. +#' from \code{/scalars/}. Default is \code{c("FD")}. #' Must match group names in the file. #' @param analysis_names Character vector. Subfolder names under #' \code{/results/} to load. Default is \code{character(0)} (none). +#' @param backend Character. Storage backend: \code{"auto"} (default), +#' \code{"hdf5"}, or \code{"tiledb"}. Auto-detection resolves to TileDB for +#' a \code{.tdb} path, or for a directory that already contains a +#' \code{scalars/} or \code{results/} subdirectory; everything else is +#' treated as HDF5. Pass \code{"tiledb"} explicitly for a TileDB store that +#' is neither. #' #' @return A \linkS4class{ModelArray} object. #' #' @seealso \linkS4class{ModelArray} for the class definition, -#' \code{\link{h5summary}} for inspecting an HDF5 file. +#' \code{\link{ModelArraySummary}} for inspecting storage. #' #' @examples #' \dontrun{ @@ -115,7 +131,8 @@ ModelArraySeed <- function(filepath, name, type = NA) { #' @export ModelArray <- function(filepath, scalar_types = c("FD"), - analysis_names = character(0)) { + analysis_names = character(0), + backend = c("auto", "hdf5", "tiledb")) { # TODO: try and use hdf5r instead of rhdf5 and delayedarray here # fn.h5 <- H5File$new(filepath, mode="a") # open; "a": creates a new file or opens an existing one for read/write @@ -127,6 +144,11 @@ ModelArray <- function(filepath, # TODO: IN THE FUTURE, THE SCALAR_TYPES AND ANALYSIS_NAMES ARE AUTOMATICALLY DETECTED # (at least detect + provide some options) + backend <- .resolve_storage_backend(filepath, backend) + if (identical(backend, "tiledb")) { + .require_tiledb_support() + } + ## scalar_data: sources <- vector("list", length(scalar_types)) scalar_data <- vector("list", length(scalar_types)) @@ -138,67 +160,11 @@ ModelArray <- function(filepath, scalar_data[[x]] <- ModelArraySeed( filepath, name = sprintf("scalars/%s/values", scalar_types[x]), - type = NA + type = NA, + backend = backend ) %>% DelayedArray::DelayedArray() - # load source filenames (column_names): prefer attribute; fallback to dataset - attrs <- rhdf5::h5readAttributes(filepath, name = sprintf("scalars/%s/values", scalar_types[x])) - colnames_attr <- attrs$column_names - if (is.null(colnames_attr)) { - # Fallback: attempt to read from dataset-based column names - # Try multiple plausible locations for compatibility across writers - paths_to_try <- c( - sprintf("scalars/%s/column_names", scalar_types[x]), - sprintf("scalars/%s/values/column_names", scalar_types[x]), - sprintf("scalars/scalars/%s/values/column_names", scalar_types[x]), - sprintf("scalars/scalars/%s/column_names", scalar_types[x]) - ) - - colnames_ds <- NULL - last_error <- NULL - for (p in paths_to_try) { - tmp <- tryCatch( - { - rhdf5::h5read(filepath, p) - }, - error = function(e) { - last_error <<- e - NULL - } - ) - if (!is.null(tmp)) { - colnames_ds <- tmp - if (grepl("^scalars/scalars/", p)) { - warning( - "Column names found at nested path '", p, "'. ", - "This is a known quirk from some converters (e.g., concifti).", - call. = FALSE - ) - } - break - } - } - if (is.null(colnames_ds)) { - stop(paste0( - "Neither attribute 'column_names' nor a dataset with column names found. Tried: ", - paste(paths_to_try, collapse = ", "), - if (!is.null(last_error)) paste0(". Last error: ", conditionMessage(last_error)) else "" - )) - } - # Ensure character vector, not list/matrix; trim potential null terminators and whitespace - if (is.list(colnames_ds)) { - colnames_ds <- unlist(colnames_ds, use.names = FALSE) - } - colnames_ds <- as.vector(colnames_ds) - colnames_ds <- as.character(colnames_ds) - # Trim any trailing NULs (hex 00) and surrounding whitespace for cross-language string compatibility - # Use escaped hex in pattern to avoid embedding a NUL in the source code - colnames_ds <- gsub("[\\x00]+$", "", colnames_ds, perl = TRUE, useBytes = TRUE) - colnames_ds <- trimws(colnames_ds) - sources[[x]] <- colnames_ds - } else { - sources[[x]] <- as.character(colnames_attr) - } + sources[[x]] <- .read_scalar_column_names(filepath, scalar_types[x], backend) # transpose scalar_data[[x]] if needed: if (dim(scalar_data[[x]])[2] == length(sources[[x]])) { @@ -230,8 +196,8 @@ ModelArray <- function(filepath, # user did not request any analyses; do not touch /results results_data <- list() } else { - # user requested analyses; check if results group exists in this .h5 file - flag_results_exist <- flagResultsGroupExistInh5(filepath) + # user requested analyses; check if results group exists in this storage + flag_results_exist <- .results_group_exists(filepath, backend) # message(flag_results_exist) if (flag_results_exist == FALSE) { results_data <- list() @@ -242,57 +208,13 @@ ModelArray <- function(filepath, for (x in seq_along(analysis_names)) { analysis_name <- analysis_names[x] - # we need to check if this subfolder exists in this .h5 file: - flag_analysis_exist <- flagAnalysisExistInh5(filepath, analysis_name = analysis_name) + # we need to check if this subfolder exists in this storage: + flag_analysis_exist <- .analysis_exists(filepath, analysis_name = analysis_name, backend = backend) if (flag_analysis_exist == FALSE) { stop(paste0("This analysis: ", analysis_name, " does not exist...")) } else { # exists - # Load column names for results: prefer attribute; fallback to dataset - attrs <- rhdf5::h5readAttributes(filepath, - name = sprintf("results/%s/results_matrix", analysis_name) - ) - names_results_matrix <- attrs$colnames - if (is.null(names_results_matrix)) { - # Fallback to dataset-based column names (similar to scalar handling) - paths_to_try <- c( - sprintf("results/%s/column_names", analysis_name), - sprintf("results/%s/results_matrix/column_names", analysis_name) - ) - colnames_ds <- NULL - last_error <- NULL - for (p in paths_to_try) { - tmp <- tryCatch( - { - rhdf5::h5read(filepath, p) - }, - error = function(e) { - last_error <<- e - NULL - } - ) - if (!is.null(tmp)) { - colnames_ds <- tmp - break - } - } - if (is.null(colnames_ds)) { - stop(paste0( - "Neither attribute 'colnames' nor a dataset with column names found for results. Tried: ", - paste(paths_to_try, collapse = ", "), - if (!is.null(last_error)) paste0(". Last error: ", conditionMessage(last_error)) else "" - )) - } - if (is.list(colnames_ds)) { - colnames_ds <- unlist(colnames_ds, use.names = FALSE) - } - colnames_ds <- as.vector(colnames_ds) - colnames_ds <- as.character(colnames_ds) - # Trim trailing NULs and whitespace - colnames_ds <- gsub("[\\x00]+$", "", colnames_ds, perl = TRUE, useBytes = TRUE) - colnames_ds <- trimws(colnames_ds) - names_results_matrix <- colnames_ds - } + names_results_matrix <- .read_result_column_names(filepath, analysis_name, backend) # names_results_matrix <- ModelArraySeed(filepath, name = sprintf( # "results/%s/has_names", analysis_name), type = NA) %>% @@ -305,7 +227,8 @@ ModelArray <- function(filepath, results_data[[x]]$results_matrix <- ModelArraySeed( filepath, name = sprintf("results/%s/results_matrix", analysis_name), - type = NA + type = NA, + backend = backend ) %>% DelayedArray::DelayedArray() if (dim(results_data[[x]]$results_matrix)[2] != length(names_results_matrix)) { @@ -322,28 +245,25 @@ ModelArray <- function(filepath, # /results//lut_col?: # LOOP OVER # OF COL OF $RESULTS_MATRIX, AND SEE IF THERE IS LUT_COL for (i_col in seq_along(names_results_matrix)) { object_name <- paste0("lut_forcol", as.character(i_col)) - flag_lut_exist <- flagObjectExistInh5( + flag_lut_exist <- .stored_object_exists( filepath, group_name = paste0("/results/", analysis_name), - object_name = object_name + object_name = object_name, + backend = backend ) if (flag_lut_exist == TRUE) { - lut <- ModelArraySeed( - filepath, - name = paste0("results/", analysis_name, "/", object_name), - type = NA - ) %>% DelayedArray::DelayedArray() + lut <- .read_result_lut(filepath, analysis_name, object_name, backend) # results_data[[x]]$lut[[i_col]] <- lut # turn values in results_matrix into factors | # HOWEVER, this also makes the entire $results_matrix into type "character".... - lut <- lut %>% as.character() - for (j_lut in seq_along(lut)) { - str_lut <- lut[j_lut] - idx_list <- results_data[[x]]$results_matrix[, i_col] %in% c(j_lut) - results_data[[x]]$results_matrix[idx_list, i_col] <- lut[j_lut] - } + lut <- as.character(lut) + results_data[[x]]$results_matrix <- .recode_result_lut_column( + results_data[[x]]$results_matrix, + i_col, + lut + ) # } else { # the lut for this column does not exist # results_data[[x]]$lut[[i_col]] <- NULL @@ -369,7 +289,8 @@ ModelArray <- function(filepath, scalars = scalar_data, results = results_data, # TODO: issue: LHS SHOULD BE THE SAME AS THE NAME IN THE H5 FILE, NOT NECESSARY CALLED "results" - path = filepath + path = filepath, + backend = backend ) } @@ -1072,18 +993,17 @@ analyseOneElement.wrap <- function(i_element, } } -#' Write outputs from element-wise statistical analysis to an HDF5 file +#' Write outputs from element-wise statistical analysis to storage #' #' @description #' Creates a group named \code{analysis_name} under \code{/results/} in the -#' HDF5 file, then writes the statistical results data.frame (i.e. for one +#' storage backend, then writes the statistical results data.frame (i.e. for one #' analysis) into it as \code{results_matrix} along with column names. #' #' @details #' The results are stored at #' \code{/results//results_matrix} with column names saved -#' as a separate dataset at -#' \code{/results//column_names}. +#' as a separate HDF5 dataset or TileDB metadata. #' #' If any column of \code{df.output} is not numeric or integer, it is #' coerced to numeric via \code{factor()} and the factor levels are saved @@ -1095,9 +1015,9 @@ analyseOneElement.wrap <- function(i_element, #' check if the message mentions "No such file or directory". Try using an #' absolute path for the \code{fn.output} argument. #' -#' @param fn.output Character. The HDF5 (\code{.h5}) filename for the output. -#' The file must already exist; use an absolute path if you encounter -#' file-not-found errors. +#' @param fn.output Character. The HDF5 (\code{.h5}) filename or TileDB +#' (\code{.tdb}) store for the output. Use an absolute path if you +#' encounter file-not-found errors. #' @param df.output A data.frame of element-wise statistical results, as #' returned by \code{\link{ModelArray.lm}}, #' \code{\link{ModelArray.gam}}, or \code{\link{ModelArray.wrap}}. @@ -1106,17 +1026,23 @@ analyseOneElement.wrap <- function(i_element, #' as the group name under \code{/results/} in the HDF5 file. #' Default is \code{"myAnalysis"}. #' @param overwrite Logical. If a group with the same \code{analysis_name} -#' already exists in the HDF5 file, whether to overwrite it (\code{TRUE}) +#' already exists in the storage backend, whether to overwrite it (\code{TRUE}) #' or skip with a warning (\code{FALSE}). Default is \code{TRUE}. +#' @param backend Character. Storage backend: \code{"auto"} (default), +#' \code{"hdf5"}, or \code{"tiledb"}. Auto-detection resolves to TileDB for +#' a \code{.tdb} path, or for a directory that already contains a +#' \code{scalars/} or \code{results/} subdirectory; everything else is +#' treated as HDF5. To create a new TileDB store at a path without a +#' \code{.tdb} suffix, pass \code{"tiledb"} explicitly. #' #' @return Invisible \code{NULL}. Called for its side effect of writing -#' results to the HDF5 file. +#' results. #' #' @seealso \code{\link{ModelArray.lm}}, \code{\link{ModelArray.gam}}, #' \code{\link{ModelArray.wrap}} which produce the \code{df.output}, #' \code{\link{results}} for reading results back from a -#' \linkS4class{ModelArray}, \code{\link{h5summary}} for inspecting what -#' has been written. +#' \linkS4class{ModelArray}, \code{\link{ModelArraySummary}} for inspecting +#' what has been written. #' #' @examples #' \dontrun{ @@ -1138,7 +1064,7 @@ analyseOneElement.wrap <- function(i_element, #' ) #' #' # Verify -#' h5summary("data.h5") +#' ModelArraySummary("data.h5") #' } #' #' @rdname writeResults @@ -1146,7 +1072,8 @@ analyseOneElement.wrap <- function(i_element, writeResults <- function(fn.output, df.output, analysis_name = "myAnalysis", - overwrite = TRUE) { + overwrite = TRUE, + backend = c("auto", "hdf5", "tiledb")) { # This is enhanced version with: 1) change to hdf5r; 2) write results with only one row for one element # check "df.output" @@ -1154,6 +1081,16 @@ writeResults <- function(fn.output, stop("Results dataset is not correct; must be data of type `data.frame`") } + backend <- .resolve_storage_backend(fn.output, backend) + if (identical(backend, "tiledb")) { + return(.write_results_tiledb( + fn.output = fn.output, + df.output = df.output, + analysis_name = analysis_name, + overwrite = overwrite + )) + } + fn.output.h5 <- hdf5r::H5File$new(fn.output, mode = "a") # open; "a": creates a new file or opens an existing one for read/write @@ -1189,37 +1126,13 @@ writeResults <- function(fn.output, results.analysis.grp <- results.grp$create_group(analysis_name) # create a subgroup called analysis_name under results.grp - # check "df.output": make sure all columns are floats (i.e. numeric) - for (i_col in seq(1, ncol(df.output), by = 1)) { - # for each column of df.output - col_class <- as.character(sapply(df.output, class)[i_col]) # class of this column - - not_numeric_or_int <- (col_class != "numeric") && - (col_class != "integer") - if (not_numeric_or_int) { - # the column class is not numeric or integer - message( - paste0( - "the column #", - as.character(i_col), - " of df.output to save: ", - "data class is not numeric or integer...fixing it" - ) - ) - - # turn into numeric && write the notes in .h5 file...: - factors <- df.output %>% - dplyr::pull(., var = i_col) %>% - factor() - df.output[, i_col] <- df.output %>% - dplyr::pull(., var = i_col) %>% - factor() %>% - as.numeric(.) # change into numeric of 1,2,3.... - - # write a LUT for this column: - results.analysis.grp[[paste0("lut_forcol", as.character(i_col))]] <- levels(factors) - # save lut to .h5/results//lut_col - } + # check "df.output": make sure all columns are floats (i.e. numeric). + # Shared with the TileDB writer so both backends derive the same LUTs. + prepared <- .prepare_results_for_storage(df.output) + df.output <- prepared$data + for (lut_name in names(prepared$luts)) { + # save lut to .h5/results//lut_forcol + results.analysis.grp[[lut_name]] <- prepared$luts[[lut_name]] } # save: diff --git a/R/ModelArray_S4Methods.R b/R/ModelArray_S4Methods.R index 109c6f1..a0f284c 100644 --- a/R/ModelArray_S4Methods.R +++ b/R/ModelArray_S4Methods.R @@ -4,8 +4,8 @@ #' @rdname ModelArray-class #' #' @description -#' Prints a summary of the ModelArray including file path, scalar dimensions, -#' and any saved analysis names. +#' Prints a summary of the ModelArray including file path, source count, +#' each scalar with its element count, and any saved analysis names. #' #' @param object A \linkS4class{ModelArray} object. #' @@ -32,21 +32,27 @@ setMethod("show", "ModelArray", function(object) { # , group_name_results="resul cat(is(object)[[1]], " located at ", path_str, "\n\n", sep = "") scalar_names <- names(scalars(object)) - for (sn in scalar_names) { - nr <- nrow(scalars(object)[[sn]]) - nc <- ncol(scalars(object)[[sn]]) - cat(format(paste0(" ", sn, ":"), justify = "left", width = 20), - nr, " elements x ", nc, " input files\n", - sep = "" - ) - } + source_counts <- lengths(sources(object)) + n_source_files <- if (length(source_counts) > 0) max(source_counts) else 0L + cat(format(" Source files:", justify = "left", width = 20), + n_source_files, "\n", + sep = "" + ) + # Element counts come from the array dimensions, which both backends carry + # in their schema — reporting them here costs no data read. + scalar_labels <- vapply(scalar_names, function(sn) { + paste0(sn, " (", nrow(scalars(object)[[sn]]), " elements)") + }, character(1)) + cat(format(" Scalars:", justify = "left", width = 20), + paste0(scalar_labels, collapse = ", "), "\n", + sep = "" + ) + analysis_names <- names(results(object)) - if (length(analysis_names) > 0) { - cat(format(" Analyses:", justify = "left", width = 20), - paste0(analysis_names, collapse = ", "), "\n", - sep = "" - ) - } + cat(format(" Analyses:", justify = "left", width = 20), + paste0(analysis_names, collapse = ", "), "\n", + sep = "" + ) }) ### Accessors for ModelArray ##### @@ -434,14 +440,16 @@ setMethod("analysisNames", "ModelArray", function(x) { #' #' @description #' Reads element metadata (e.g., greyordinates for CIFTI data, or fixel/voxel -#' coordinate information) from the HDF5 file if present. The function searches -#' for known metadata dataset names (\code{"greyordinates"}, \code{"fixels"}, -#' \code{"voxels"}) at the top level of the HDF5 file. +#' coordinate information) from the HDF5 file or TileDB store if present. The +#' function searches for known metadata dataset names (\code{"greyordinates"}, +#' \code{"fixels"}, \code{"voxels"}) at the top level of the storage backend. +#' It uses the backend retained when the \code{ModelArray} was constructed, so +#' an explicit constructor choice is preserved. #' #' @param x A \linkS4class{ModelArray} object. #' #' @return A matrix or data.frame of element metadata if found, or \code{NULL} -#' if no known metadata dataset exists in the HDF5 file. +#' if no known metadata dataset exists. #' #' @seealso \code{\link{nElements}}, \code{\link{scalars}} #' @@ -464,6 +472,32 @@ setMethod("elementMetadata", "ModelArray", function(x) { # Try known metadata dataset names metadata_paths <- c("greyordinates", "fixels", "voxels") + backend <- .modelarray_backend(x, filepath) + if (identical(backend, "tiledb")) { + .require_tiledb_support() + for (p in metadata_paths) { + result <- tryCatch( + { + uri <- .modelarray_path(filepath, p) + if (!.tiledb_object_exists(uri, "ARRAY")) { + NULL + } else { + result <- as.matrix(DelayedArray::DelayedArray( + TileDBArray::TileDBArraySeed(uri, attr = "values") + )) + dimnames(result) <- NULL + result + } + }, + error = function(e) NULL + ) + if (!is.null(result)) { + return(result) + } + } + return(NULL) + } + for (p in metadata_paths) { result <- tryCatch( rhdf5::h5read(filepath, p), diff --git a/R/analyse-helpers.R b/R/analyse-helpers.R index e8186d5..1b1db68 100644 --- a/R/analyse-helpers.R +++ b/R/analyse-helpers.R @@ -399,6 +399,132 @@ # This replaces duplicated logic across analyseOneElement.lm, # analyseOneElement.gam, and analyseOneElement.wrap. +#' Return the storage path associated with one scalar +#' @noRd +.scalar_storage_path <- function(modelarray, scalar) { + paths <- modelarray@path + if (length(paths) == 0L) { + return(NA_character_) + } + if (!is.null(names(paths)) && scalar %in% names(paths)) { + return(paths[[scalar]]) + } + paths[[1]] +} + + +#' Return TRUE when any attached scalar is TileDB-backed +#' @noRd +.has_tiledb_attached_scalar <- function(ctx) { + any(vapply(ctx$attached_scalars, function(sname) { + path <- .scalar_storage_path(ctx$modelarray, sname) + if (is.na(path) || is.null(path)) { + return(FALSE) + } + identical(.modelarray_backend(ctx$modelarray, path), "tiledb") + }, logical(1))) +} + + +#' Resolve the scalar row cache block size for a context +#' +#' TileDB row-at-a-time reads are expensive in element-wise analyses. When a +#' TileDB-backed scalar participates in the model, materialize rows in chunks +#' before dispatching to workers. Set option `ModelArray.tiledb_read_block_size` +#' to override the memory-bounded default, or set it to 0 to disable. Set +#' `ModelArray.tiledb_read_block_mb` to tune the default memory target. +#' @noRd +.scalar_read_block_size <- function(ctx) { + if (!.has_tiledb_attached_scalar(ctx)) { + return(NULL) + } + + block_size <- getOption("ModelArray.tiledb_read_block_size", NULL) + if (is.null(block_size)) { + block_mb <- getOption("ModelArray.tiledb_read_block_mb", 512) + if (!is.numeric(block_mb) || length(block_mb) != 1L || is.na(block_mb) || block_mb <= 0) { + stop("Option ModelArray.tiledb_read_block_mb must be a positive single number") + } + n_cols <- max(vapply(ctx$attached_scalars, function(sname) { + ncol(scalars(ctx$modelarray)[[sname]]) + }, numeric(1)), 1L) + bytes_per_row <- max(1L, n_cols) * max(1L, length(ctx$attached_scalars)) * 8L + target_bytes <- block_mb * 1024^2 + block_size <- max(1L, min(.Machine$integer.max, floor(target_bytes / bytes_per_row))) + } + + if (!is.numeric(block_size) || length(block_size) != 1L || is.na(block_size)) { + stop("Option ModelArray.tiledb_read_block_size must be a single number") + } + block_size <- as.integer(min(block_size, .Machine$integer.max)) + if (block_size <= 0L) { + return(NULL) + } + max(block_size, 1L) +} + + +#' Add a materialized scalar row cache to an analysis context +#' +#' The cache is intentionally chunk-local. In forked parallel execution this +#' lets the parent process perform the TileDB I/O once and workers reuse the +#' in-memory block through copy-on-write. +#' @noRd +.with_scalar_row_cache <- function(ctx, element.subset) { + if (is.null(.scalar_read_block_size(ctx))) { + return(ctx) + } + + element.subset <- as.integer(element.subset) + if (length(element.subset) == 0L) { + return(ctx) + } + + cache <- stats::setNames( + vector("list", length(ctx$attached_scalars)), + ctx$attached_scalars + ) + for (sname in ctx$attached_scalars) { + cache[[sname]] <- as.matrix( + scalars(ctx$modelarray)[[sname]][element.subset, , drop = FALSE] + ) + } + + ctx$scalar_row_cache <- cache + ctx$scalar_row_cache_elements <- element.subset + ctx$scalar_row_cache_contiguous <- length(element.subset) == 1L || + all(diff(element.subset) == 1L) + ctx$scalar_row_cache_first <- element.subset[[1]] + ctx$scalar_row_cache_position <- stats::setNames( + seq_along(element.subset), + as.character(element.subset) + ) + ctx +} + + +#' Return a cached or storage-backed scalar row +#' @noRd +.read_scalar_row <- function(ctx, sname, i_element) { + cache <- ctx$scalar_row_cache + if (!is.null(cache) && !is.null(cache[[sname]])) { + if (isTRUE(ctx$scalar_row_cache_contiguous)) { + pos <- as.integer(i_element - ctx$scalar_row_cache_first + 1L) + } else { + # Single-bracket: `[[` on a named vector raises "subscript out of + # bounds" for an absent name rather than returning NULL, which would + # defeat the fall-through to a direct storage read below. + pos <- unname(ctx$scalar_row_cache_position[as.character(i_element)]) + } + if (!is.na(pos) && pos >= 1L && pos <= nrow(cache[[sname]])) { + return(as.vector(cache[[sname]][pos, ])) + } + } + + as.vector(scalars(ctx$modelarray)[[sname]][i_element, ]) +} + + #' Assemble per-element data.frame from precomputed context #' #' Reads scalar rows from the ModelArray, applies precomputed reorder @@ -423,7 +549,7 @@ #' @noRd .assemble_element_data <- function(i_element, ctx, num.subj.lthr) { # Read the response scalar row — the only mandatory per-element I/O - response_vals <- scalars(ctx$modelarray)[[ctx$scalar]][i_element, ] + response_vals <- .read_scalar_row(ctx, ctx$scalar, i_element) # Start the validity mask with the response scalar valid_mask <- is.finite(response_vals) @@ -434,7 +560,7 @@ other_scalars <- setdiff(ctx$attached_scalars, ctx$scalar) for (sname in other_scalars) { - s_vals <- scalars(ctx$modelarray)[[sname]][i_element, ] + s_vals <- .read_scalar_row(ctx, sname, i_element) # Apply precomputed reorder index (or use as-is if NULL) reorder_idx <- ctx$predictor_reorder[[sname]] @@ -593,6 +719,23 @@ fits } + +#' Print progress context for chunked element dispatch +#' @noRd +.report_iteration_progress <- function(iteration, total_iterations, + chunk_start, chunk_end, + total_elements, pbar) { + if (!isTRUE(pbar) || total_iterations <= 1L) { + return(invisible(NULL)) + } + message( + "iteration ", iteration, "/", total_iterations, + " (elements ", chunk_start, "-", chunk_end, + " of ", total_elements, ")" + ) + invisible(NULL) +} + # P-value correction ---- #' Correct p-values for a set of terms and append corrected columns #' @@ -622,6 +765,274 @@ } # Streaming writes ---- +#' Validate a gzip compression level used by a streaming writer +#' @noRd +.validate_stream_compression_level <- function(compression_level, argument) { + if (!is.numeric(compression_level) || + length(compression_level) != 1L || + is.na(compression_level) || + !is.finite(compression_level) || + compression_level != floor(compression_level) || + compression_level < 0L || + compression_level > 9L) { + stop(argument, " must be a single integer between 0 and 9", call. = FALSE) + } + as.integer(compression_level) +} + + +#' Create a dense TileDB values array for streaming +#' @noRd +.create_tiledb_values_array <- function(uri, + n_rows, + n_cols, + chunk_rows, + storage_mode, + compression_level) { + .require_tiledb_support() + + supported_modes <- c("logical", "integer", "double") + if (!is.character(storage_mode) || + length(storage_mode) != 1L || + !storage_mode %in% supported_modes) { + stop( + "TileDB streaming supports storage modes: ", + paste(supported_modes, collapse = ", "), + call. = FALSE + ) + } + + parent_uri <- dirname(uri) + if (dir.exists(parent_uri) || .tiledb_object_exists(parent_uri)) { + unlink(parent_uri, recursive = TRUE, force = TRUE) + } + dir.create(parent_uri, recursive = TRUE, showWarnings = FALSE) + + dimensions <- list( + tiledb::tiledb_dim( + name = "d1", + domain = c(1L, as.integer(n_rows)), + tile = as.integer(chunk_rows), + type = "INT32" + ), + tiledb::tiledb_dim( + name = "d2", + domain = c(1L, as.integer(n_cols)), + tile = as.integer(n_cols), + type = "INT32" + ) + ) + domain <- tiledb::tiledb_domain(dims = dimensions) + + gzip_filter <- tiledb::tiledb_filter("GZIP") + gzip_filter <- tiledb::tiledb_filter_set_option( + gzip_filter, + "COMPRESSION_LEVEL", + compression_level + ) + value_filters <- tiledb::tiledb_filter_list(list(gzip_filter)) + value_attribute <- tiledb::tiledb_attr( + name = "values", + type = tiledb::r_to_tiledb_type(vector(storage_mode)), + filter_list = value_filters + ) + schema <- tiledb::tiledb_array_schema( + domain = domain, + attrs = list(value_attribute), + cell_order = "COL_MAJOR", + tile_order = "COL_MAJOR", + sparse = FALSE + ) + tiledb::tiledb_array_create(uri, schema) + + .write_tiledb_array_metadata(uri, list(type = storage_mode)) + invisible(NULL) +} + + +#' Write one matrix block to a dense TileDB values array +#' @noRd +.write_tiledb_values_block <- function(uri, row_idx, n_cols, block) { + array <- tiledb::tiledb_array(uri) + on.exit(try(tiledb::tiledb_array_close(array), silent = TRUE), add = TRUE) + query <- tiledb::tiledb_query(array, "WRITE") + query <- tiledb::tiledb_query_set_subarray( + query, + as.integer(c(range(row_idx), 1L, n_cols)), + "INT32" + ) + query <- tiledb::tiledb_query_set_layout(query, "COL_MAJOR") + query <- tiledb::tiledb_query_set_buffer( + query, + "values", + as.vector(block) + ) + query <- tiledb::tiledb_query_submit(query) + query <- tiledb::tiledb_query_finalize(query) + if (!identical(tiledb::tiledb_query_status(query), "COMPLETE")) { + stop("TileDB block write did not complete", call. = FALSE) + } + invisible(NULL) +} + + +#' Initialize an incremental writer for /scalars datasets +#' @noRd +.init_scalar_stream_writer <- function(write_scalar_name, + write_scalar_file, + n_rows, + column_names, + flush_every = 1000L, + storage_mode = "double", + compression_level = 4L) { + compression_level <- .validate_stream_compression_level( + compression_level, + "write_scalar_compression_level" + ) + if (!is.character(write_scalar_file) || + length(write_scalar_file) != 1L || + is.na(write_scalar_file) || + !nzchar(write_scalar_file)) { + stop("write_scalar_file must be a non-empty character path") + } + backend <- .resolve_storage_backend(write_scalar_file, "auto") + n_cols <- length(column_names) + chunk_rows <- min(flush_every, n_rows) + + if (identical(backend, "tiledb")) { + scalar_uri <- .storage_child_path( + .modelarray_path(write_scalar_file, "scalars"), + write_scalar_name, + "write_scalar_name" + ) + values_uri <- .modelarray_path(scalar_uri, "values") + .create_tiledb_values_array( + uri = values_uri, + n_rows = n_rows, + n_cols = n_cols, + chunk_rows = chunk_rows, + storage_mode = storage_mode, + compression_level = compression_level + ) + return(list( + backend = "tiledb", + file = write_scalar_file, + scalar_name = write_scalar_name, + values_uri = values_uri, + column_names = as.character(column_names), + storage_mode = storage_mode, + n_cols = n_cols, + write_row_cursor = 1L + )) + } + + if (!file.exists(write_scalar_file)) { + rhdf5::h5createFile(write_scalar_file) + } + scalar_group <- paste0("scalars/", write_scalar_name) + h5_write <- hdf5r::H5File$new(write_scalar_file, mode = "a") + if (!h5_write$exists("scalars")) { + h5_write$create_group("scalars") + } + scalars_group <- h5_write$open("scalars") + if (scalars_group$exists(write_scalar_name)) { + scalars_group$link_delete(write_scalar_name) + } + scalars_group$create_group(write_scalar_name) + h5_write$close_all() + + dataset_path <- paste0(scalar_group, "/values") + rhdf5::h5createDataset( + file = write_scalar_file, + dataset = dataset_path, + dims = c(n_rows, n_cols), + storage.mode = storage_mode, + chunk = c(chunk_rows, n_cols), + level = compression_level + ) + + list( + backend = "hdf5", + file = write_scalar_file, + scalar_name = write_scalar_name, + dataset_path = dataset_path, + column_names = as.character(column_names), + storage_mode = storage_mode, + n_cols = n_cols, + write_row_cursor = 1L + ) +} + + +#' Append one block to an incremental /scalars writer +#' @noRd +.scalar_stream_write_block <- function(writer, block) { + if (is.null(writer)) { + return(writer) + } + block <- as.matrix(block) + if (ncol(block) != writer$n_cols) { + stop("Scalar block does not match the initialized writer column count") + } + storage.mode(block) <- writer$storage_mode + row_idx <- writer$write_row_cursor:(writer$write_row_cursor + nrow(block) - 1L) + + if (identical(writer$backend, "tiledb")) { + .write_tiledb_values_block( + uri = writer$values_uri, + row_idx = row_idx, + n_cols = writer$n_cols, + block = block + ) + } else { + rhdf5::h5write( + obj = block, + file = writer$file, + name = writer$dataset_path, + index = list(row_idx, seq_len(writer$n_cols)) + ) + } + + writer$write_row_cursor <- writer$write_row_cursor + nrow(block) + writer +} + + +#' Finalize an incremental /scalars writer +#' @noRd +.finalize_scalar_stream_writer <- function(writer) { + if (is.null(writer)) { + return(invisible(NULL)) + } + if (identical(writer$backend, "tiledb")) { + .write_tiledb_array_metadata( + writer$values_uri, + list( + column_names = jsonlite::toJSON( + writer$column_names, + auto_unbox = TRUE + ) + ) + ) + return(invisible(NULL)) + } + + .write_hdf5_attribute( + writer$file, + writer$dataset_path, + "column_names", + writer$column_names + ) + rhdf5::h5write( + obj = writer$column_names, + file = writer$file, + name = paste0("scalars/", writer$scalar_name, "/column_names") + ) + rhdf5::h5closeAll() + invisible(NULL) +} + + #' Initialize an incremental writer for /results datasets #' @noRd .init_results_stream_writer <- function(write_results_name, @@ -634,16 +1045,64 @@ if (is.null(write_results_name)) { return(NULL) } - if (!is.character(write_results_name) || length(write_results_name) != 1L || write_results_name == "") { + if (!is.character(write_results_name) || + length(write_results_name) != 1L || + is.na(write_results_name) || + !nzchar(write_results_name)) { stop("write_results_name must be a non-empty character string when provided") } - if (is.null(write_results_file) || !is.character(write_results_file) || length(write_results_file) != 1L) { + if (is.null(write_results_file) || + !is.character(write_results_file) || + length(write_results_file) != 1L || + is.na(write_results_file) || + !nzchar(write_results_file)) { stop("write_results_file must be a single character path when write_results_name is provided") } if (!is.numeric(flush_every) || length(flush_every) != 1L || flush_every <= 0) { stop("write_results_flush_every must be a positive integer") } flush_every <- as.integer(flush_every) + compression_level <- .validate_stream_compression_level( + compression_level, + "write_results_compression_level" + ) + + backend <- .resolve_storage_backend(write_results_file, "auto") + if (identical(backend, "tiledb")) { + .require_tiledb_support() + analysis_uri <- .storage_child_path( + .modelarray_path(write_results_file, "results"), + write_results_name, + "write_results_name" + ) + values_uri <- .modelarray_path(analysis_uri, "results_matrix") + .create_tiledb_values_array( + uri = values_uri, + n_rows = n_rows, + n_cols = length(column_names), + chunk_rows = min(flush_every, n_rows), + storage_mode = storage_mode, + compression_level = compression_level + ) + return(list( + backend = "tiledb", + file = write_results_file, + analysis_name = write_results_name, + values_uri = values_uri, + column_names = as.character(column_names), + storage_mode = storage_mode, + n_cols = length(column_names), + write_row_cursor = 1L + )) + } + + if (!identical(backend, "hdf5")) { + stop( + "Unsupported write_results_file backend: ", + backend, + call. = FALSE + ) + } if (!file.exists(write_results_file)) { rhdf5::h5createFile(write_results_file) @@ -669,11 +1128,13 @@ dims = c(n_rows, length(column_names)), storage.mode = storage_mode, chunk = c(chunk_rows, length(column_names)), - level = as.integer(compression_level) + level = compression_level ) list( + backend = "hdf5", file = write_results_file, + analysis_name = write_results_name, dataset_path = dataset_path, names_path = paste0("results/", write_results_name, "/column_names"), column_names = as.character(column_names), @@ -689,8 +1150,23 @@ if (is.null(writer)) { return(writer) } + if (ncol(block_df) != writer$n_cols) { + stop("Result block does not match the initialized writer column count") + } block <- as.matrix(block_df) row_idx <- writer$write_row_cursor:(writer$write_row_cursor + nrow(block) - 1L) + if (identical(writer$backend, "tiledb")) { + storage.mode(block) <- writer$storage_mode + .write_tiledb_values_block( + uri = writer$values_uri, + row_idx = row_idx, + n_cols = writer$n_cols, + block = block + ) + writer$write_row_cursor <- writer$write_row_cursor + nrow(block) + return(writer) + } + rhdf5::h5write( obj = block, file = writer$file, @@ -708,6 +1184,19 @@ if (is.null(writer)) { return(invisible(NULL)) } + if (identical(writer$backend, "tiledb")) { + .write_tiledb_array_metadata( + writer$values_uri, + list( + column_names = jsonlite::toJSON( + writer$column_names, + auto_unbox = TRUE + ) + ) + ) + return(invisible(NULL)) + } + rhdf5::h5write( obj = writer$column_names, file = writer$file, diff --git a/R/analyse.R b/R/analyse.R index 09e3620..ce095f1 100644 --- a/R/analyse.R +++ b/R/analyse.R @@ -75,11 +75,12 @@ #' skips. Default: \code{"stop"}. #' @param write_results_name Optional character. If provided, results are #' incrementally written to -#' \code{results//results_matrix} in the HDF5 file +#' \code{results//results_matrix} in the output #' specified by \code{write_results_file}. -#' @param write_results_file Optional character. HDF5 file path for -#' incremental result writes. Required when \code{write_results_name} -#' is provided. +#' @param write_results_file Optional character. HDF5 file path or TileDB +#' store path for incremental result writes. Required when +#' \code{write_results_name} is provided. The backend is detected from the +#' path, so use a \code{.tdb} suffix to stream into a TileDB store. #' @param write_results_flush_every Positive integer. Number of elements #' per write block. Default 1000. #' @param write_results_storage_mode Character. Storage mode for HDF5 @@ -88,7 +89,7 @@ #' level for HDF5 writes. Default 4. #' @param return_output Logical. If \code{TRUE} (default), return the #' combined data.frame. If \code{FALSE}, return \code{invisible(NULL)}; -#' useful when writing large outputs directly to HDF5. +#' useful when writing outputs directly to storage. #' @param ... Additional arguments passed to \code{\link[stats]{lm}}. #' #' @return A tibble with one row per element. The first column is @@ -271,16 +272,26 @@ ModelArray.lm <- function(formula, data, phenotypes, scalar = NULL, element.subs # Start model looping ---- fits_all <- if (need_full_df) vector("list", length(element.subset)) else NULL chunk_size <- if (is.null(writer)) length(element.subset) else as.integer(write_results_flush_every) + read_block_size <- .scalar_read_block_size(ctx) + if (!is.null(read_block_size)) { + chunk_size <- min(chunk_size, read_block_size) + } chunk_starts <- seq(1L, length(element.subset), by = chunk_size) - for (chunk_start in chunk_starts) { + for (chunk_idx in seq_along(chunk_starts)) { + chunk_start <- chunk_starts[[chunk_idx]] chunk_end <- min(chunk_start + chunk_size - 1L, length(element.subset)) chunk_elements <- element.subset[chunk_start:chunk_end] + chunk_ctx <- .with_scalar_row_cache(ctx, chunk_elements) + .report_iteration_progress( + chunk_idx, length(chunk_starts), + chunk_start, chunk_end, length(element.subset), pbar + ) ## Main loop now passes ctx instead of formula/modelarray/phenotypes/scalar chunk_fits <- .parallel_dispatch( chunk_elements, analyseOneElement.lm, n_cores, pbar, - ctx = ctx, + ctx = chunk_ctx, var.terms = var.terms, var.model = var.model, num.subj.lthr = num.subj.lthr, num.stat.output = num.stat.output, flag_initiate = FALSE, on_error = on_error, @@ -310,13 +321,14 @@ ModelArray.lm <- function(formula, data, phenotypes, scalar = NULL, element.subs df_out <- .correct_pvalues(df_out, list.terms, correct.p.value.terms, var.terms) df_out <- .correct_pvalues(df_out, "model", correct.p.value.model, var.model) - # Rewrite corrected results to HDF5 if needed ---- + # Rewrite corrected results if streaming was used ---- if (!is.null(writer) && (need_term_correction || need_model_correction)) { writeResults( fn.output = write_results_file, df.output = df_out, analysis_name = write_results_name, - overwrite = TRUE + overwrite = TRUE, + backend = writer$backend ) } @@ -666,16 +678,26 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar = NULL, element.sub # Start model looping ---- fits_all <- if (need_full_df) vector("list", length(element.subset)) else NULL chunk_size <- if (is.null(writer)) length(element.subset) else as.integer(write_results_flush_every) + read_block_size <- .scalar_read_block_size(ctx) + if (!is.null(read_block_size)) { + chunk_size <- min(chunk_size, read_block_size) + } chunk_starts <- seq(1L, length(element.subset), by = chunk_size) - for (chunk_start in chunk_starts) { + for (chunk_idx in seq_along(chunk_starts)) { + chunk_start <- chunk_starts[[chunk_idx]] chunk_end <- min(chunk_start + chunk_size - 1L, length(element.subset)) chunk_elements <- element.subset[chunk_start:chunk_end] + chunk_ctx <- .with_scalar_row_cache(ctx, chunk_elements) + .report_iteration_progress( + chunk_idx, length(chunk_starts), + chunk_start, chunk_end, length(element.subset), pbar + ) # Main loop passes ctx chunk_fits <- .parallel_dispatch( chunk_elements, analyseOneElement.gam, n_cores, pbar, - ctx = ctx, + ctx = chunk_ctx, var.smoothTerms = var.smoothTerms, var.parametricTerms = var.parametricTerms, var.model = var.model, @@ -751,16 +773,39 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar = NULL, element.sub reduced.model.column_names <- reduced.model.outputs_initiator$column_names reduced.model.num.stat.output <- length(reduced.model.column_names) - reduced.model.fits <- .parallel_dispatch( - element.subset, analyseOneElement.gam, n_cores, pbar, - ctx = reduced_ctx, - var.smoothTerms = c(), var.parametricTerms = c(), - var.model = c("adj.r.squared"), - num.subj.lthr = num.subj.lthr, - num.stat.output = reduced.model.num.stat.output, - flag_initiate = FALSE, flag_sse = TRUE, - ... - ) + reduced.read_block_size <- .scalar_read_block_size(reduced_ctx) + reduced.chunk_size <- if (is.null(reduced.read_block_size)) { + length(element.subset) + } else { + min(length(element.subset), reduced.read_block_size) + } + reduced.chunk_starts <- seq(1L, length(element.subset), by = reduced.chunk_size) + reduced.model.fits <- vector("list", length(element.subset)) + + for (reduced.chunk_idx in seq_along(reduced.chunk_starts)) { + reduced.chunk_start <- reduced.chunk_starts[[reduced.chunk_idx]] + reduced.chunk_end <- min( + reduced.chunk_start + reduced.chunk_size - 1L, + length(element.subset) + ) + reduced.chunk_elements <- element.subset[reduced.chunk_start:reduced.chunk_end] + reduced.chunk_ctx <- .with_scalar_row_cache(reduced_ctx, reduced.chunk_elements) + .report_iteration_progress( + reduced.chunk_idx, length(reduced.chunk_starts), + reduced.chunk_start, reduced.chunk_end, length(element.subset), pbar + ) + + reduced.model.fits[reduced.chunk_start:reduced.chunk_end] <- .parallel_dispatch( + reduced.chunk_elements, analyseOneElement.gam, n_cores, pbar, + ctx = reduced.chunk_ctx, + var.smoothTerms = c(), var.parametricTerms = c(), + var.model = c("adj.r.squared"), + num.subj.lthr = num.subj.lthr, + num.stat.output = reduced.model.num.stat.output, + flag_initiate = FALSE, flag_sse = TRUE, + ... + ) + } result_mat_reduced <- do.call(rbind, reduced.model.fits) col_list_reduced <- lapply( @@ -797,7 +842,8 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar = NULL, element.sub fn.output = write_results_file, df.output = df_out, analysis_name = write_results_name, - overwrite = TRUE + overwrite = TRUE, + backend = writer$backend ) } @@ -843,10 +889,12 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar = NULL, element.sub #' vector of results for that element. #' @param write_scalar_name Optional character. If provided, selected #' output columns are written into -#' \code{scalars//values} in the HDF5 file specified -#' by \code{write_scalar_file}. -#' @param write_scalar_file Optional character. HDF5 output file path. -#' Required when \code{write_scalar_name} is provided. +#' \code{scalars//values} in the output specified by +#' \code{write_scalar_file}. +#' @param write_scalar_file Optional character. HDF5 file path or TileDB store +#' path. Required when \code{write_scalar_name} is provided. The backend is +#' detected from the path, so use a \code{.tdb} suffix to write a TileDB +#' store. #' @param write_scalar_columns Optional character or integer vector #' selecting which output columns to save as scalar values. If #' \code{NULL} (default), uses all output columns except @@ -857,9 +905,10 @@ ModelArray.gam <- function(formula, data, phenotypes, scalar = NULL, element.sub #' @param write_scalar_flush_every Positive integer. Elements per write #' block for scalar writes. Default 1000. #' @param write_scalar_storage_mode Character. Storage mode for scalar -#' writes (e.g. \code{"double"}). Default \code{"double"}. +#' writes (e.g. \code{"double"}). TileDB supports \code{"logical"}, +#' \code{"integer"}, and \code{"double"}. Default \code{"double"}. #' @param write_scalar_compression_level Integer 0--9. Gzip compression -#' level for scalar writes. Default 4. +#' level for HDF5 and TileDB scalar writes. Default 4. #' @param ... Additional arguments forwarded to \code{FUN}. #' #' @return If \code{flag_initiate = TRUE}, a list with one component: @@ -978,6 +1027,7 @@ ModelArray.wrap <- function(FUN, data, phenotypes, scalar, element.subset = NULL ) writing_scalar <- !is.null(write_scalar_name) + writer_scalar <- NULL if (writing_scalar) { if (!is.character(write_scalar_name) || length(write_scalar_name) != 1L || write_scalar_name == "") { stop("write_scalar_name must be a non-empty character string when provided") @@ -985,7 +1035,9 @@ ModelArray.wrap <- function(FUN, data, phenotypes, scalar, element.subset = NULL if (is.null(write_scalar_file) || !is.character(write_scalar_file) || length(write_scalar_file) != 1L) { stop("write_scalar_file must be a single character path when write_scalar_name is provided") } - if (!is.numeric(write_scalar_flush_every) || length(write_scalar_flush_every) != 1L || write_scalar_flush_every <= 0) { + if (!is.numeric(write_scalar_flush_every) || + length(write_scalar_flush_every) != 1L || + write_scalar_flush_every <= 0) { stop("write_scalar_flush_every must be a positive integer") } write_scalar_flush_every <- as.integer(write_scalar_flush_every) @@ -1011,38 +1063,19 @@ ModelArray.wrap <- function(FUN, data, phenotypes, scalar, element.subset = NULL stop("length(write_scalar_column_names) must equal number of selected write_scalar_columns") } - # Initialize scalar output dataset - if (!file.exists(write_scalar_file)) { - rhdf5::h5createFile(write_scalar_file) - } - scalar_grp <- paste0("scalars/", write_scalar_name) - h5_write <- hdf5r::H5File$new(write_scalar_file, mode = "a") - if (!h5_write$exists("scalars")) { - h5_write$create_group("scalars") - } - scalars_grp_obj <- h5_write$open("scalars") - if (scalars_grp_obj$exists(write_scalar_name)) { - scalars_grp_obj$link_delete(write_scalar_name) - } - scalars_grp_obj$create_group(write_scalar_name) - h5_write$close_all() - - dataset_path <- paste0(scalar_grp, "/values") - chunk_rows <- min(write_scalar_flush_every, length(element.subset)) - rhdf5::h5createDataset( - file = write_scalar_file, - dataset = dataset_path, - dims = c(length(element.subset), length(scalar_col_idx)), - storage.mode = write_scalar_storage_mode, - chunk = c(chunk_rows, length(scalar_col_idx)), - level = as.integer(write_scalar_compression_level) + writer_scalar <- .init_scalar_stream_writer( + write_scalar_name = write_scalar_name, + write_scalar_file = write_scalar_file, + n_rows = length(element.subset), + column_names = write_scalar_column_names, + flush_every = write_scalar_flush_every, + storage_mode = write_scalar_storage_mode, + compression_level = write_scalar_compression_level ) } # Start model looping ---- fits_all <- if (return_output) vector("list", length(element.subset)) else NULL - write_row_cursor <- 1L - chunk_size <- length(element.subset) if (writing_scalar) { chunk_size <- min(chunk_size, write_scalar_flush_every) @@ -1050,17 +1083,27 @@ ModelArray.wrap <- function(FUN, data, phenotypes, scalar, element.subset = NULL if (!is.null(writer_results)) { chunk_size <- min(chunk_size, as.integer(write_results_flush_every)) } + read_block_size <- .scalar_read_block_size(ctx) + if (!is.null(read_block_size)) { + chunk_size <- min(chunk_size, read_block_size) + } chunk_starts <- seq(1L, length(element.subset), by = chunk_size) - for (chunk_start in chunk_starts) { + for (chunk_idx in seq_along(chunk_starts)) { + chunk_start <- chunk_starts[[chunk_idx]] chunk_end <- min(chunk_start + chunk_size - 1L, length(element.subset)) chunk_elements <- element.subset[chunk_start:chunk_end] + chunk_ctx <- .with_scalar_row_cache(ctx, chunk_elements) + .report_iteration_progress( + chunk_idx, length(chunk_starts), + chunk_start, chunk_end, length(element.subset), pbar + ) ## Main loop passes ctx ---- chunk_fits <- .parallel_dispatch( chunk_elements, analyseOneElement.wrap, n_cores, pbar, user_fun = FUN, - ctx = ctx, + ctx = chunk_ctx, num.subj.lthr = num.subj.lthr, num.stat.output = num.stat.output, flag_initiate = FALSE, on_error = on_error, @@ -1082,32 +1125,12 @@ ModelArray.wrap <- function(FUN, data, phenotypes, scalar, element.subset = NULL if (writing_scalar) { block <- as.matrix(chunk_df[, scalar_col_idx, drop = FALSE]) - row_idx <- write_row_cursor:(write_row_cursor + nrow(block) - 1L) - rhdf5::h5write( - obj = block, - file = write_scalar_file, - name = dataset_path, - index = list(row_idx, seq_len(ncol(block))) - ) - write_row_cursor <- write_row_cursor + nrow(block) + writer_scalar <- .scalar_stream_write_block(writer_scalar, block) } } - # Finalize H5 writing ---- - if (writing_scalar) { - rhdf5::h5writeAttribute( - attr = write_scalar_column_names, - h5obj = write_scalar_file, - name = "column_names", - h5loc = dataset_path - ) - rhdf5::h5write( - obj = write_scalar_column_names, - file = write_scalar_file, - name = paste0("scalars/", write_scalar_name, "/column_names") - ) - rhdf5::h5closeAll() - } + # Finalize streaming writes ---- + .finalize_scalar_stream_writer(writer_scalar) .finalize_results_stream_writer(writer_results) if (!return_output) { diff --git a/R/merge.R b/R/merge.R index 9638ac8..ec768e3 100644 --- a/R/merge.R +++ b/R/merge.R @@ -160,6 +160,7 @@ mergeModelArrays <- function(modelarrays, phenotypes_list, merge_on) { combined_scalars <- list() combined_sources <- list() combined_paths <- character() + combined_backends <- character() for (i in seq_along(modelarrays)) { ma <- modelarrays[[i]] @@ -193,6 +194,7 @@ mergeModelArrays <- function(modelarrays, phenotypes_list, merge_on) { colnames(combined_scalars[[sn]]) <- unified_source_id combined_sources[[sn]] <- unified_source_id combined_paths[sn] <- ma@path[1] + combined_backends[sn] <- .modelarray_backend(ma, ma@path[1]) } } @@ -253,7 +255,8 @@ mergeModelArrays <- function(modelarrays, phenotypes_list, merge_on) { scalars = combined_scalars, sources = combined_sources, results = list(), - path = combined_paths + path = combined_paths, + backend = combined_backends ) list(data = combined_ma, phenotypes = merged_phen) diff --git a/R/utils.R b/R/utils.R index 93a9acd..462ee6d 100644 --- a/R/utils.R +++ b/R/utils.R @@ -33,6 +33,590 @@ flagAnalysisExistInh5 <- function(fn_h5, analysis_name) { } +#' Write an HDF5 attribute with hdf5r's current API +#' @noRd +.write_hdf5_attribute <- function(filepath, object_path, attr_name, value) { + h5 <- hdf5r::H5File$new(filepath, mode = "a") + on.exit(h5$close_all(), add = TRUE) + object <- h5[[object_path]] + # create_attr() errors if the attribute is already there, unlike the + # rhdf5::h5writeAttribute() it replaces, which overwrote silently. + if (object$attr_exists(attr_name)) { + object$attr_delete(attr_name) + } + object$create_attr(attr_name, robj = value) + invisible(NULL) +} + + +#' Does a directory look like a ModelArray storage tree? +#' +#' Filesystem-only on purpose: this runs for HDF5 paths too, so it must not +#' reach for the optional TileDB packages. +#' @noRd +.looks_like_modelarray_store <- function(filepath) { + dir.exists(file.path(filepath, "scalars")) || + dir.exists(file.path(filepath, "results")) +} + + +#' Resolve a storage backend from a path and user preference +#' @noRd +.resolve_storage_backend <- function(filepath, backend = c("auto", "hdf5", "tiledb")) { + backend <- match.arg(backend) + if (backend != "auto") { + return(backend) + } + + # A `.tdb` suffix declares intent, whether or not the store exists yet, so it + # is also how a *new* TileDB store gets created without passing `backend`. + if (grepl("\\.tdb/?$", filepath, ignore.case = TRUE)) { + return("tiledb") + } + # Otherwise only claim a directory when it actually holds a store. Writers + # delete the group they are about to replace, so a mistyped output directory + # must not be mistaken for one. + if (dir.exists(filepath) && .looks_like_modelarray_store(filepath)) { + return("tiledb") + } + "hdf5" +} + + +#' Return the backend retained by a ModelArray object +#' @noRd +.modelarray_backend <- function(x, filepath = x@path[1]) { + stored_backend <- tryCatch( + methods::slot(x, "backend"), + error = function(e) character(0) + ) + + if (length(stored_backend) > 1L && length(x@path) == length(stored_backend)) { + path_index <- match(filepath, x@path) + if (!is.na(path_index)) { + stored_backend <- stored_backend[path_index] + } + } + + # unname(): mergeModelArrays() stores one backend per scalar name, and every + # caller compares the result with identical(., "tiledb") — a surviving name + # makes that FALSE, which would silently disable TileDB handling on any + # merged object. + stored_backend <- unname(stored_backend[1]) + if ( + length(stored_backend) == 0L || + is.na(stored_backend) || + !stored_backend %in% c("hdf5", "tiledb") + ) { + return(.resolve_storage_backend(filepath, "auto")) + } + stored_backend +} + + +#' Require TileDB packages for TileDB-backed operations +#' @noRd +.require_tiledb_support <- function() { + missing <- c( + if (!requireNamespace("tiledb", quietly = TRUE)) "tiledb", + if (!requireNamespace("TileDBArray", quietly = TRUE)) "TileDBArray", + if (!requireNamespace("jsonlite", quietly = TRUE)) "jsonlite" + ) + if (length(missing) > 0) { + stop( + "TileDB support requires the optional package(s): ", + paste(missing, collapse = ", "), + ". Install them before using a TileDB-backed ModelArray.", + call. = FALSE + ) + } +} + + +#' Join path segments for local paths and TileDB URIs +#' @noRd +.modelarray_path <- function(...) { + parts <- c(...) + parts <- gsub("/+$", "", parts) + parts <- gsub("^/+", "", parts) + first <- c(...)[[1]] + first <- gsub("/+$", "", first) + paste(c(first, parts[-1]), collapse = "/") +} + + +#' Validate a storage object name as one safe path component +#' @noRd +.validate_storage_component <- function(value, argument) { + if ( + !is.character(value) || + length(value) != 1L || + is.na(value) || + !nzchar(value) + ) { + stop(argument, " must be a non-empty character string", call. = FALSE) + } + if (value %in% c(".", "..") || grepl("[/\\\\]", value)) { + stop( + argument, + " must be a single storage name, cannot be '.' or '..', and cannot contain '/' or '\\'", + call. = FALSE + ) + } + value +} + + +#' Build a storage child path and verify local-path containment +#' @noRd +.storage_child_path <- function(parent, component, argument) { + component <- .validate_storage_component(component, argument) + child <- .modelarray_path(parent, component) + + is_uri <- grepl("^[[:alpha:]][[:alnum:]+.-]*://", parent) + if (!is_uri) { + normalized_parent <- normalizePath(parent, winslash = "/", mustWork = FALSE) + normalized_child <- if (file.exists(child) || dir.exists(child)) { + normalizePath(child, winslash = "/", mustWork = FALSE) + } else { + file.path(normalized_parent, component) + } + parent_prefix <- paste0(sub("/+$", "", normalized_parent), "/") + if (!startsWith(normalized_child, parent_prefix)) { + stop(argument, " resolves outside its storage parent", call. = FALSE) + } + } + + child +} + + +#' Normalize vectors read from metadata or companion datasets +#' @noRd +.clean_name_vector <- function(x) { + if (is.null(x)) { + return(NULL) + } + if (is.list(x) && !is.data.frame(x)) { + x <- unlist(x, use.names = FALSE) + } + x <- as.vector(x) + x <- as.character(x) + x <- gsub("[\\x00]+$", "", x, perl = TRUE, useBytes = TRUE) + trimws(x) +} + + +#' Read names from HDF5 attributes, falling back to datasets +#' @noRd +.read_hdf5_names <- function(filepath, attr_path, attr_name, dataset_paths, label) { + attrs <- rhdf5::h5readAttributes(filepath, name = attr_path) + names_attr <- attrs[[attr_name]] + if (!is.null(names_attr)) { + return(.clean_name_vector(names_attr)) + } + + colnames_ds <- NULL + last_error <- NULL + for (p in dataset_paths) { + tmp <- tryCatch( + rhdf5::h5read(filepath, p), + error = function(e) { + last_error <<- e + NULL + } + ) + if (!is.null(tmp)) { + colnames_ds <- tmp + if (grepl("^scalars/scalars/", p)) { + warning( + "Column names found at nested path '", p, "'. ", + "This is a known quirk from some converters (e.g., concifti).", + call. = FALSE + ) + } + break + } + } + + if (is.null(colnames_ds)) { + stop(paste0( + "Neither attribute '", attr_name, "' nor a dataset with ", label, " found. Tried: ", + paste(dataset_paths, collapse = ", "), + if (!is.null(last_error)) paste0(". Last error: ", conditionMessage(last_error)) else "" + )) + } + .clean_name_vector(colnames_ds) +} + + +#' Return a TileDB object type, or INVALID when lookup fails +#' @noRd +.tiledb_object_type <- function(uri) { + .require_tiledb_support() + tryCatch(tiledb::tiledb_object_type(uri), error = function(e) "INVALID") +} + + +#' Does a TileDB object exist? +#' @noRd +.tiledb_object_exists <- function(uri, type = NULL) { + object_type <- .tiledb_object_type(uri) + if (is.null(type)) { + return(toupper(object_type) %in% c("ARRAY", "GROUP")) + } + identical(toupper(object_type), toupper(type)) +} + + +#' Convert TileDB metadata values to a character vector +#' @noRd +.parse_tiledb_names <- function(x) { + if (is.null(x)) { + return(NULL) + } + if (is.raw(x)) { + x <- rawToChar(x) + } + if (is.list(x) && length(x) == 1L) { + x <- x[[1]] + } + x <- .clean_name_vector(x) + if (length(x) == 1L && jsonlite::validate(x)) { + x <- jsonlite::fromJSON(x) + } + .clean_name_vector(x) +} + + +#' Read metadata from a TileDB array +#' @noRd +.read_tiledb_array_metadata <- function(uri, key) { + .require_tiledb_support() + if (!.tiledb_object_exists(uri, "ARRAY")) { + return(NULL) + } + arr <- tiledb::tiledb_array(uri, query_type = "READ", keep_open = TRUE) + on.exit(try(tiledb::tiledb_array_close(arr), silent = TRUE), add = TRUE) + tryCatch(tiledb::tiledb_get_metadata(arr, key), error = function(e) NULL) +} + + +#' Read metadata from a TileDB group +#' @noRd +.read_tiledb_group_metadata <- function(uri, key) { + .require_tiledb_support() + if (!.tiledb_object_exists(uri, "GROUP")) { + return(NULL) + } + grp <- tiledb::tiledb_group(uri, type = "READ") + on.exit(try(tiledb::tiledb_group_close(grp), silent = TRUE), add = TRUE) + tryCatch(tiledb::tiledb_group_get_metadata(grp, key), error = function(e) NULL) +} + + +#' Read the "values" attribute from a 1-D TileDB array +#' @noRd +.read_tiledb_values_array <- function(uri) { + .require_tiledb_support() + if (!.tiledb_object_exists(uri, "ARRAY")) { + return(NULL) + } + arr <- tiledb::tiledb_array( + uri, + query_type = "READ", + attrs = "values", + return_as = "data.frame" + ) + on.exit(try(tiledb::tiledb_array_close(arr), silent = TRUE), add = TRUE) + values <- tryCatch(arr[], error = function(e) NULL) + if (is.null(values)) { + return(NULL) + } + if (is.data.frame(values) && "values" %in% names(values)) { + return(.clean_name_vector(values[["values"]])) + } + if (is.list(values) && "values" %in% names(values)) { + return(.clean_name_vector(values[["values"]])) + } + .clean_name_vector(values) +} + + +#' Return dimensions for a dense TileDB array +#' @noRd +.tiledb_array_dims <- function(uri) { + .require_tiledb_support() + arr <- tiledb::tiledb_array(uri, query_type = "READ", keep_open = TRUE) + on.exit(try(tiledb::tiledb_array_close(arr), silent = TRUE), add = TRUE) + + dom <- tiledb::domain(tiledb::schema(arr)) + ndim <- tiledb::tiledb_ndim(dom) + vapply(seq_len(ndim), function(i) { + extent <- tiledb::tiledb_array_get_non_empty_domain_from_index(arr, i) + as.integer(as.numeric(extent[2]) - as.numeric(extent[1]) + 1) + }, integer(1)) +} + + +#' Read scalar column names for HDF5 or TileDB backends +#' @noRd +.read_scalar_column_names <- function(filepath, scalar_type, backend) { + if (identical(backend, "hdf5")) { + return(.read_hdf5_names( + filepath = filepath, + attr_path = sprintf("scalars/%s/values", scalar_type), + attr_name = "column_names", + dataset_paths = c( + sprintf("scalars/%s/column_names", scalar_type), + sprintf("scalars/%s/values/column_names", scalar_type), + sprintf("scalars/scalars/%s/values/column_names", scalar_type), + sprintf("scalars/scalars/%s/column_names", scalar_type) + ), + label = "column names" + )) + } + + .require_tiledb_support() + value_uri <- .modelarray_path(filepath, "scalars", scalar_type, "values") + scalar_uri <- .modelarray_path(filepath, "scalars", scalar_type) + names <- .parse_tiledb_names(.read_tiledb_array_metadata(value_uri, "column_names")) + if (is.null(names)) { + names <- .parse_tiledb_names(.read_tiledb_group_metadata(scalar_uri, "column_names")) + } + if (is.null(names)) { + names <- .read_tiledb_values_array(.modelarray_path(scalar_uri, "column_names")) + } + if (is.null(names)) { + stop( + "No TileDB column names found for scalar '", scalar_type, "'. ", + "Expected column_names metadata on ", value_uri, + " or a values array at ", .modelarray_path(scalar_uri, "column_names"), + call. = FALSE + ) + } + names +} + + +#' Read result column names for HDF5 or TileDB backends +#' @noRd +.read_result_column_names <- function(filepath, analysis_name, backend) { + if (identical(backend, "hdf5")) { + return(.read_hdf5_names( + filepath = filepath, + attr_path = sprintf("results/%s/results_matrix", analysis_name), + attr_name = "colnames", + dataset_paths = c( + sprintf("results/%s/column_names", analysis_name), + sprintf("results/%s/results_matrix/column_names", analysis_name) + ), + label = "column names for results" + )) + } + + .require_tiledb_support() + matrix_uri <- .modelarray_path(filepath, "results", analysis_name, "results_matrix") + analysis_uri <- .modelarray_path(filepath, "results", analysis_name) + names <- .parse_tiledb_names(.read_tiledb_array_metadata(matrix_uri, "colnames")) + if (is.null(names)) { + names <- .parse_tiledb_names(.read_tiledb_array_metadata(matrix_uri, "column_names")) + } + if (is.null(names)) { + names <- .parse_tiledb_names(.read_tiledb_group_metadata(analysis_uri, "column_names")) + } + if (is.null(names)) { + names <- .read_tiledb_values_array(.modelarray_path(analysis_uri, "column_names")) + } + if (is.null(names)) { + stop("No TileDB column names found for analysis '", analysis_name, "'.", call. = FALSE) + } + names +} + + +#' Test for result groups/arrays in either supported backend +#' @noRd +.results_group_exists <- function(filepath, backend) { + if (identical(backend, "hdf5")) { + return(flagResultsGroupExistInh5(filepath)) + } + .tiledb_object_exists(.modelarray_path(filepath, "results"), "GROUP") || + dir.exists(.modelarray_path(filepath, "results")) +} + + +#' Test for a named analysis in either supported backend +#' @noRd +.analysis_exists <- function(filepath, analysis_name, backend) { + if (identical(backend, "hdf5")) { + return(flagAnalysisExistInh5(filepath, analysis_name = analysis_name)) + } + analysis_uri <- .modelarray_path(filepath, "results", analysis_name) + matrix_uri <- .modelarray_path(analysis_uri, "results_matrix") + .tiledb_object_exists(analysis_uri) || + dir.exists(analysis_uri) || + .tiledb_object_exists(matrix_uri, "ARRAY") +} + + +#' Test for a named object in either supported backend +#' @noRd +.stored_object_exists <- function(filepath, group_name, object_name, backend) { + if (identical(backend, "hdf5")) { + return(flagObjectExistInh5(filepath, group_name = group_name, object_name = object_name)) + } + group_uri <- .modelarray_path(filepath, gsub("^/", "", group_name)) + object_uri <- .modelarray_path(group_uri, object_name) + .tiledb_object_exists(object_uri, "ARRAY") || + !is.null(.read_tiledb_array_metadata(.modelarray_path(group_uri, "results_matrix"), object_name)) || + !is.null(.read_tiledb_group_metadata(group_uri, object_name)) +} + + +#' Read a result look-up table for either supported backend +#' @noRd +.read_result_lut <- function(filepath, analysis_name, object_name, backend) { + if (identical(backend, "hdf5")) { + if (!flagObjectExistInh5( + filepath, + group_name = paste0("/results/", analysis_name), + object_name = object_name + )) { + return(NULL) + } + return(ModelArraySeed( + filepath, + name = paste0("results/", analysis_name, "/", object_name), + type = NA, + backend = "hdf5" + ) %>% + DelayedArray::DelayedArray() %>% + as.character()) + } + + analysis_uri <- .modelarray_path(filepath, "results", analysis_name) + matrix_uri <- .modelarray_path(analysis_uri, "results_matrix") + lut <- .parse_tiledb_names(.read_tiledb_array_metadata(matrix_uri, object_name)) + if (is.null(lut)) { + lut <- .parse_tiledb_names(.read_tiledb_group_metadata(analysis_uri, object_name)) + } + if (is.null(lut)) { + lut <- .read_tiledb_values_array(.modelarray_path(analysis_uri, object_name)) + } + lut +} + + +#' Recode one results matrix column from stored LUT indices to labels +#' @noRd +.recode_result_lut_column <- function(results_matrix, i_col, lut) { + lut <- as.character(lut) + stored_values <- as.character(as.vector(results_matrix[, i_col])) + lut_index <- match(stored_values, as.character(seq_along(lut))) + recoded_values <- stored_values + recoded_values[!is.na(lut_index)] <- lut[lut_index[!is.na(lut_index)]] + results_matrix[, i_col] <- recoded_values + results_matrix +} + + +#' Coerce result data.frame columns to storage-friendly numeric columns +#' @noRd +.prepare_results_for_storage <- function(df.output) { + luts <- list() + for (i_col in seq_len(ncol(df.output))) { + col <- df.output[[i_col]] + not_numeric_or_int <- !is.numeric(col) && !is.integer(col) + if (not_numeric_or_int) { + message( + paste0( + "the column #", + as.character(i_col), + " of df.output to save: ", + "data class is not numeric or integer...fixing it" + ) + ) + + factors <- factor(col) + df.output[[i_col]] <- as.numeric(factors) + luts[[paste0("lut_forcol", as.character(i_col))]] <- levels(factors) + } + } + list(data = df.output, luts = luts) +} + + +#' Write metadata to a TileDB array, ignoring close failures +#' @noRd +.write_tiledb_array_metadata <- function(uri, values) { + .require_tiledb_support() + arr <- tiledb::tiledb_array(uri, query_type = "WRITE", keep_open = TRUE) + on.exit(try(tiledb::tiledb_array_close(arr), silent = TRUE), add = TRUE) + for (key in names(values)) { + value <- values[[key]] + if (inherits(value, "json")) { + value <- as.character(value) + } + tiledb::tiledb_put_metadata(arr, key, value) + } + invisible(NULL) +} + + +#' Write ModelArray results to a TileDB store +#' @noRd +.write_results_tiledb <- function(fn.output, df.output, analysis_name, overwrite) { + .require_tiledb_support() + + if ( + !is.character(fn.output) || + length(fn.output) != 1L || + is.na(fn.output) || + !nzchar(fn.output) + ) { + stop("fn.output must be a non-empty character path", call. = FALSE) + } + results_uri <- .modelarray_path(fn.output, "results") + analysis_uri <- .storage_child_path(results_uri, analysis_name, "analysis_name") + + prepared <- .prepare_results_for_storage(df.output) + df.output <- prepared$data + + if (!dir.exists(fn.output)) { + dir.create(fn.output, recursive = TRUE, showWarnings = FALSE) + } + + matrix_uri <- .modelarray_path(analysis_uri, "results_matrix") + analysis_exists <- dir.exists(analysis_uri) || .tiledb_object_exists(analysis_uri) + + if (analysis_exists && !overwrite) { + warning(paste0(analysis_name, " exists but not to overwrite!")) + return(invisible(NULL)) + } + if (analysis_exists && overwrite) { + unlink(analysis_uri, recursive = TRUE, force = TRUE) + } + + dir.create(analysis_uri, recursive = TRUE, showWarnings = FALSE) + TileDBArray::writeTileDBArray( + as.matrix(df.output), + path = matrix_uri, + attr = "values" + ) + + metadata <- c( + list(column_names = jsonlite::toJSON(as.character(colnames(df.output)), auto_unbox = TRUE)), + stats::setNames( + lapply(prepared$luts, jsonlite::toJSON, auto_unbox = TRUE), + names(prepared$luts) + ) + ) + .write_tiledb_array_metadata(matrix_uri, metadata) + invisible(NULL) +} + + #' print the additional arguments settings #' @param FUN The function, e.g. mgcv::gam, without "()" #' @param argu_name The argument name of the function @@ -99,8 +683,9 @@ check_validity_correctPValue <- function(correct.list, name.correct.list, "p.value was not included in ", name.var.list, ", so not to perform its p.value corrections" - ) - ) # TODO: why this warning comes out after ModelArray.aModel is done? + ), + immediate. = TRUE + ) } } } @@ -527,29 +1112,40 @@ bind_cols_check_emptyTibble <- function(a, b) { } -#' Summarize an HDF5 file without loading a full ModelArray +#' Summarize ModelArray storage without loading a full ModelArray #' #' @description -#' Reads the HDF5 file structure and returns a summary of available scalars, +#' Reads the storage structure and returns a summary of available scalars, #' their dimensions, and any saved analyses. Useful for inspecting large files #' without constructing a full \linkS4class{ModelArray} object. #' #' @details -#' This function opens the HDF5 file read-only via \code{\link[rhdf5]{h5ls}}, -#' inspects the group structure under \code{/scalars/} and \code{/results/}, -#' and closes the file. It does not load any data into memory. The returned +#' For HDF5, this function opens the file read-only via +#' \code{\link[rhdf5]{h5ls}}. For TileDB, it inspects arrays under +#' \code{/scalars/} and \code{/results/}. It does not load full data matrices +#' into memory. The returned #' object has a \code{print} method that displays a formatted summary. #' -#' @param filepath Character. Path to an HDF5 (\code{.h5}) file. +#' \code{h5summary()} is a backward-compatible alias for +#' \code{ModelArraySummary()}. Despite the historical name, it supports both +#' HDF5 files and TileDB stores. #' -#' @return An object of class \code{"h5summary"}, which is a list with -#' components: +#' @param filepath Character. Path to an HDF5 (\code{.h5}) file or TileDB +#' (\code{.tdb}) store. +#' @param backend Character. Storage backend: \code{"auto"} (default), +#' \code{"hdf5"}, or \code{"tiledb"}. +#' +#' @return \code{ModelArraySummary()} returns an object of class +#' \code{c("ModelArraySummary", "h5summary")}. \code{h5summary()} returns +#' an object of class \code{"h5summary"} for backward compatibility. Both +#' are lists with components: #' \describe{ #' \item{scalars}{A data.frame with columns \code{name}, #' \code{nElements}, and \code{nInputFiles}.} #' \item{analyses}{Character vector of analysis names found under #' \code{/results/}.} #' \item{filepath}{The input filepath.} +#' \item{backend}{The resolved storage backend.} #' } #' #' @seealso \code{\link{ModelArray}} for loading the full object, @@ -557,20 +1153,120 @@ bind_cols_check_emptyTibble <- function(a, b) { #' #' @examples #' \dontrun{ -#' h5summary("path/to/data.h5") +#' ModelArraySummary("path/to/data.h5") +#' ModelArraySummary("path/to/store.tdb") #' #' # Inspect before deciding which scalars to load -#' info <- h5summary("path/to/data.h5") +#' info <- ModelArraySummary("path/to/data.h5") #' info$scalars$name #' ma <- ModelArray("path/to/data.h5", scalar_types = info$scalars$name) +#' +#' # Historical alias, still supported +#' h5summary("path/to/data.h5") #' } #' -#' @rdname h5summary +#' @rdname ModelArraySummary +#' @export +ModelArraySummary <- function(filepath, backend = c("auto", "hdf5", "tiledb")) { + structure( + .modelarray_summary_impl(filepath, backend = backend), + class = c("ModelArraySummary", "h5summary") + ) +} + +#' @rdname ModelArraySummary #' @export -h5summary <- function(filepath) { +h5summary <- function(filepath, backend = c("auto", "hdf5", "tiledb")) { + structure( + .modelarray_summary_impl(filepath, backend = backend), + class = "h5summary" + ) +} + +#' @noRd +.modelarray_summary_impl <- function(filepath, backend = c("auto", "hdf5", "tiledb")) { if (!file.exists(filepath)) { stop("File not found: ", filepath) } + backend <- .resolve_storage_backend(filepath, backend) + if (identical(backend, "tiledb")) { + .require_tiledb_support() + scalars_dir <- .modelarray_path(filepath, "scalars") + scalar_candidates <- if (dir.exists(scalars_dir)) { + basename(list.dirs(scalars_dir, recursive = FALSE, full.names = TRUE)) + } else { + character(0) + } + scalar_names <- scalar_candidates[vapply( + scalar_candidates, + function(scalar_name) { + .tiledb_object_exists( + .modelarray_path(filepath, "scalars", scalar_name, "values"), + "ARRAY" + ) + }, + logical(1) + )] + scalar_info <- do.call(rbind, lapply(scalar_names, function(scalar_name) { + values_uri <- .modelarray_path(filepath, "scalars", scalar_name, "values") + dims <- .tiledb_array_dims(values_uri) + source_names <- .read_scalar_column_names(filepath, scalar_name, backend) + n_input_files <- length(source_names) + if (dims[2] == n_input_files) { + n_elements <- dims[1] + } else if (dims[1] == n_input_files) { + n_elements <- dims[2] + } else { + stop( + "TileDB scalar '", scalar_name, "' dimensions ", + paste(dims, collapse = " x "), + " do not match its ", n_input_files, " stored column names", + call. = FALSE + ) + } + data.frame( + name = scalar_name, + nElements = as.integer(n_elements), + nInputFiles = as.integer(n_input_files), + stringsAsFactors = FALSE + ) + })) + if (is.null(scalar_info)) { + scalar_info <- data.frame( + name = character(0), + nElements = integer(0), + nInputFiles = integer(0), + stringsAsFactors = FALSE + ) + } + + results_dir <- .modelarray_path(filepath, "results") + analysis_candidates <- if (dir.exists(results_dir)) { + basename(list.dirs(results_dir, recursive = FALSE, full.names = TRUE)) + } else { + character(0) + } + analyses <- analysis_candidates[vapply( + analysis_candidates, + function(analysis_name) { + .tiledb_object_exists( + .modelarray_path(filepath, "results", analysis_name, "results_matrix"), + "ARRAY" + ) + }, + logical(1) + )] + + return( + list( + scalars = scalar_info, + analyses = analyses, + filepath = filepath, + backend = backend + ) + ) + } + listing <- rhdf5::h5ls(filepath) rhdf5::h5closeAll() @@ -601,29 +1297,40 @@ h5summary <- function(filepath) { ] analyses <- result_groups$name - structure( - list( - scalars = scalar_info, - analyses = analyses, - filepath = filepath - ), - class = "h5summary" + list( + scalars = scalar_info, + analyses = analyses, + filepath = filepath, + backend = backend ) } -#' @rdname h5summary +#' @rdname ModelArraySummary #' -#' @param x An \code{h5summary} object as returned by -#' \code{\link{h5summary}}. +#' @param x A \code{ModelArraySummary} or \code{h5summary} object. #' @param ... Additional arguments (currently ignored). #' #' @return Invisible \code{x}. Called for its side effect of printing a #' human-readable summary to the console. #' +#' @method print ModelArraySummary +#' @export +print.ModelArraySummary <- function(x, ...) { + .print_modelarray_summary(x, ...) +} + +#' @rdname ModelArraySummary #' @method print h5summary #' @export print.h5summary <- function(x, ...) { - cat("H5 file:", x$filepath, "\n\n") + .print_modelarray_summary(x, ...) +} + +#' @noRd +.print_modelarray_summary <- function(x, ...) { + backend <- if (is.null(x$backend)) "hdf5" else x$backend + label <- if (identical(backend, "tiledb")) "TileDB store" else "H5 file" + cat(label, ": ", x$filepath, "\n\n", sep = "") if (nrow(x$scalars) > 0) { cat("Scalars:\n") for (i in seq_len(nrow(x$scalars))) { diff --git a/README.Rmd b/README.Rmd index 49cb379..403c9bc 100644 --- a/README.Rmd +++ b/README.Rmd @@ -46,7 +46,7 @@ That's it. ## Why ModelArray? -**It scales.** ModelArray uses [HDF5](https://www.hdfgroup.org/solutions/hdf5/) for on-disk storage and [DelayedArray](https://bioconductor.org/packages/DelayedArray/) for lazy access. +**It scales.** ModelArray uses [HDF5](https://www.hdfgroup.org/solutions/hdf5/) by default, and can also read TileDB stores when the optional [TileDB](https://tiledb.com/) R packages are installed. Both backends use [DelayedArray](https://bioconductor.org/packages/DelayedArray/) for lazy access. An ABCC dMRI dataset has ~350,000 voxels × ~26,000-sessions - about 291 GB if loaded in its entirety - never enters RAM. ModelArray reads one element at a time, so memory usage stays flat regardless of dataset size. @@ -55,7 +55,7 @@ ModelArray reads one element at a time, so memory usage stays flat regardless of - **GAMs** (`ModelArray.gam()`) with penalized splines for capturing nonlinear effects, such as lifespan trajectories that accelerate or decelerate across development - **Any R modeling function** (`ModelArray.wrap()`) — if you can write a function that takes a data frame and returns a row of statistics, ModelArray will run it across your entire dataset -**It works with what you have.** ModelArray is modality-agnostic — the same code works for fixel-wise data, voxel-wise data, and surface-based greyordinate data. The companion tool [ModelArrayIO](https://github.com/PennLINC/ModelArrayIO) handles conversion from `.mif`, NIfTI, and CIFTI formats into the HDF5 file that ModelArray expects. +**It works with what you have.** ModelArray is modality-agnostic — the same code works for fixel-wise data, voxel-wise data, and surface-based greyordinate data. The companion tool [ModelArrayIO](https://github.com/PennLINC/ModelArrayIO) handles conversion from `.mif`, NIfTI, and CIFTI formats into HDF5 files or TileDB stores that ModelArray can read. **It's been peer-reviewed.** ModelArray was published in *NeuroImage* in 2023 and has been used in studies of lifespan brain development. @@ -63,7 +63,7 @@ ModelArray reads one element at a time, so memory usage stays flat regardless of
-![Overview](man/figures/overview_structure.svg) +[![Overview](https://raw.githubusercontent.com/PennLINC/ModelArrayIO/main/docs/_static/overview_structure.png)](https://github.com/PennLINC/ModelArrayIO/blob/main/docs/_static/overview_structure.png)
@@ -81,6 +81,17 @@ The short version: devtools::install_github("PennLINC/ModelArray") ``` +Optional TileDB support: + +```r +install.packages(c("jsonlite", "tiledb")) + +if (!requireNamespace("BiocManager", quietly = TRUE)) { + install.packages("BiocManager") +} +BiocManager::install("TileDBArray") +``` + Prefer a container? A [Docker/Singularity image](https://hub.docker.com/r/pennlinc/modelarray_confixel) with ModelArray and ModelArrayIO pre-installed is available — useful for HPC clusters where you can't install system libraries. ## Documentation diff --git a/README.md b/README.md index d5c44d8..30d9d41 100644 --- a/README.md +++ b/README.md @@ -42,7 +42,9 @@ That’s it. `results` is a data frame with estimates, *t*-statistics, ## Why ModelArray? **It scales.** ModelArray uses -[HDF5](https://www.hdfgroup.org/solutions/hdf5/) for on-disk storage and +[HDF5](https://www.hdfgroup.org/solutions/hdf5/) by default, and can +also read TileDB stores when the optional [TileDB](https://tiledb.com/) +R packages are installed. Both backends use [DelayedArray](https://bioconductor.org/packages/DelayedArray/) for lazy access. An ABCC dMRI dataset has ~350,000 voxels × ~26,000-sessions - about 291 GB if loaded in its entirety - never enters RAM. ModelArray @@ -62,8 +64,8 @@ dataset size. same code works for fixel-wise data, voxel-wise data, and surface-based greyordinate data. The companion tool [ModelArrayIO](https://github.com/PennLINC/ModelArrayIO) handles -conversion from `.mif`, NIfTI, and CIFTI formats into the HDF5 file that -ModelArray expects. +conversion from `.mif`, NIfTI, and CIFTI formats into HDF5 files or +TileDB stores that ModelArray can read. **It’s been peer-reviewed.** ModelArray was published in *NeuroImage* in 2023 and has been used in studies of lifespan brain development. @@ -72,10 +74,7 @@ ModelArray expects.
-
-Overview - -
+[![Overview](https://raw.githubusercontent.com/PennLINC/ModelArrayIO/main/docs/_static/overview_structure.png)](https://github.com/PennLINC/ModelArrayIO/blob/main/docs/_static/overview_structure.png)
@@ -97,6 +96,17 @@ The short version: devtools::install_github("PennLINC/ModelArray") ``` +Optional TileDB support: + +``` r +install.packages(c("jsonlite", "tiledb")) + +if (!requireNamespace("BiocManager", quietly = TRUE)) { + install.packages("BiocManager") +} +BiocManager::install("TileDBArray") +``` + Prefer a container? A [Docker/Singularity image](https://hub.docker.com/r/pennlinc/modelarray_confixel) with ModelArray and ModelArrayIO pre-installed is available — useful for HPC diff --git a/_pkgdown.yml b/_pkgdown.yml index 399a2fb..c1ad1ee 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -21,10 +21,10 @@ navbar: href: articles/elements.html - text: "Splitting Element IDs" href: articles/element-splitting.html - - text: "Understanding the HDF5 File" - href: articles/hdf5-format.html - - text: "HDF5 for Large-Scale Analyses" - href: articles/hdf5-large-analyses.html + - text: "Understanding ModelArray Storage" + href: articles/modelarray-storage.html + - text: "Large-Scale Analyses" + href: articles/large-scale-analyses.html case-studies: text: Case Studies menu: @@ -33,7 +33,7 @@ navbar: - text: "Modelling in Depth" href: articles/modelling.html - text: "Exploring Data with Convenience Functions" - href: articles/exploring-h5.html + href: articles/exploring-modelarray-data.html - text: "-------" - text: Resources - text: "Containers" @@ -85,6 +85,7 @@ reference: - title: "Utilities" contents: - writeResults + - ModelArraySummary - h5summary - numElementsTotal - gen_gamFormula_fxSmooth @@ -95,5 +96,8 @@ redirects: - ["articles/wrap_function.html", "articles/modelling.html"] - ["articles/voxel-wise_data.html", "articles/modelling.html"] - ["articles/basic_r_intro.html", "articles/installations.html"] - - ["articles/faq.html", "articles/exploring-h5.html"] - - ["articles/debugging.html", "articles/exploring-h5.html"] \ No newline at end of file + - ["articles/hdf5-format.html", "articles/modelarray-storage.html"] + - ["articles/hdf5-large-analyses.html", "articles/large-scale-analyses.html"] + - ["articles/exploring-h5.html", "articles/exploring-modelarray-data.html"] + - ["articles/faq.html", "articles/exploring-modelarray-data.html"] + - ["articles/debugging.html", "articles/exploring-modelarray-data.html"] diff --git a/man/ModelArray-class.Rd b/man/ModelArray-class.Rd index df25d8d..3cd63fe 100644 --- a/man/ModelArray-class.Rd +++ b/man/ModelArray-class.Rd @@ -7,7 +7,12 @@ \alias{show,ModelArray-method} \title{ModelArray class} \usage{ -ModelArray(filepath, scalar_types = c("FD"), analysis_names = character(0)) +ModelArray( + filepath, + scalar_types = c("FD"), + analysis_names = character(0), + backend = c("auto", "hdf5", "tiledb") +) \S4method{show}{ModelArray}(object) } @@ -22,19 +27,19 @@ to the console. A ModelArray wraps one or more element-wise scalar matrices (e.g., FD, FC, log_FC for fixel data) read lazily via \pkg{DelayedArray}, along with any previously saved analysis results. The object holds references to the -underlying HDF5 file and reads data on demand, making it suitable for +underlying storage and reads data on demand, making it suitable for large-scale neuroimaging datasets. -Prints a summary of the ModelArray including file path, scalar dimensions, -and any saved analysis names. +Prints a summary of the ModelArray including file path, source count, +each scalar with its element count, and any saved analysis names. } \details{ ModelArray is an S4 class that represents element-wise scalar data and -associated statistical results backed by an HDF5 file on disk. +associated statistical results backed by an HDF5 file or TileDB store on disk. -Each scalar in the HDF5 file is stored at \code{/scalars//values} +Each scalar is stored at \code{/scalars//values} as a matrix of elements (rows) by source files (columns). Source filenames -are read from HDF5 attributes or companion datasets. Analysis results, if +are read from storage metadata or companion datasets. Analysis results, if present, live under \code{/results//results_matrix}. ModelArray objects are typically created with the \code{\link{ModelArray}} @@ -55,7 +60,10 @@ Each matrix has elements as rows and source files as columns.} list containing at minimum \code{results_matrix} (a \link[DelayedArray:DelayedArray-class]{DelayedArray::DelayedArray}).} -\item{\code{path}}{Character. Path(s) to the HDF5 file(s) on disk.} +\item{\code{path}}{Character. Path(s) to the HDF5 file(s) or TileDB store(s) on disk.} + +\item{\code{backend}}{Character. Resolved storage backend(s) for \code{path}: +\code{"hdf5"} or \code{"tiledb"}.} }} \seealso{ diff --git a/man/ModelArray-package.Rd b/man/ModelArray-package.Rd index cd0cd64..ca7aa85 100644 --- a/man/ModelArray-package.Rd +++ b/man/ModelArray-package.Rd @@ -7,12 +7,13 @@ \description{ The ModelArray package provides an S4 class and associated methods for performing massively univariate statistical analyses on element-wise -(fixel, voxel, or vertex) neuroimaging data stored in HDF5 files. +(fixel, voxel, or vertex) neuroimaging data stored in HDF5 files or +TileDB stores. } \details{ The core workflow is: \enumerate{ -\item Inspect an HDF5 file with \code{\link{h5summary}()} +\item Inspect storage with \code{\link{ModelArraySummary}()} \item Load data with \code{\link{ModelArray}()} \item Fit models with \code{\link{ModelArray.lm}}, \code{\link{ModelArray.gam}}, or \code{\link{ModelArray.wrap}} @@ -40,7 +41,7 @@ analysis of fixel-wise data. \emph{NeuroImage}, 271, 120037. \linkS4class{ModelArray}, \code{\link{ModelArray}}, \code{\link{ModelArray.lm}}, \code{\link{ModelArray.gam}}, \code{\link{ModelArray.wrap}}, \code{\link{mergeModelArrays}}, -\code{\link{h5summary}} +\code{\link{ModelArraySummary}} } \author{ \strong{Maintainer}: Matthew Cieslak \email{Matthew.Cieslak@pennmedicine.upenn.edu} diff --git a/man/ModelArray.Rd b/man/ModelArray.Rd index 72c2563..2f30288 100644 --- a/man/ModelArray.Rd +++ b/man/ModelArray.Rd @@ -2,33 +2,45 @@ % Please edit documentation in R/ModelArray_Constructor.R \name{ModelArray} \alias{ModelArray} -\title{Load element-wise data from an HDF5 file} +\title{Load element-wise data from an HDF5 file or TileDB store} \usage{ -ModelArray(filepath, scalar_types = c("FD"), analysis_names = character(0)) +ModelArray( + filepath, + scalar_types = c("FD"), + analysis_names = character(0), + backend = c("auto", "hdf5", "tiledb") +) } \arguments{ \item{filepath}{Character. Path to an existing HDF5 (\code{.h5}) -file containing element-wise scalar data.} +file or TileDB (\code{.tdb}) store containing element-wise scalar data.} \item{scalar_types}{Character vector. Names of scalar groups to read -from \code{/scalars/} in the HDF5 file. Default is \code{c("FD")}. +from \code{/scalars/}. Default is \code{c("FD")}. Must match group names in the file.} \item{analysis_names}{Character vector. Subfolder names under \code{/results/} to load. Default is \code{character(0)} (none).} + +\item{backend}{Character. Storage backend: \code{"auto"} (default), +\code{"hdf5"}, or \code{"tiledb"}. Auto-detection resolves to TileDB for +a \code{.tdb} path, or for a directory that already contains a +\code{scalars/} or \code{results/} subdirectory; everything else is +treated as HDF5. Pass \code{"tiledb"} explicitly for a TileDB store that +is neither.} } \value{ A \linkS4class{ModelArray} object. } \description{ Reads scalar matrices and (optionally) saved analysis results from -an HDF5 file and returns a \linkS4class{ModelArray} object. +an HDF5 file or TileDB store and returns a \linkS4class{ModelArray} object. } \details{ The constructor reads each scalar listed in \code{scalar_types} from \code{/scalars//values}, wrapping them as \link[DelayedArray:DelayedArray-class]{DelayedArray::DelayedArray} objects. Source filenames are extracted -from HDF5 attributes or companion datasets. +from storage metadata or companion datasets. If \code{analysis_names} is non-empty, saved results are loaded from \code{/results//results_matrix}. @@ -47,5 +59,5 @@ ma } \seealso{ \linkS4class{ModelArray} for the class definition, -\code{\link{h5summary}} for inspecting an HDF5 file. +\code{\link{ModelArraySummary}} for inspecting storage. } diff --git a/man/ModelArray.gam.Rd b/man/ModelArray.gam.Rd index e08e336..072a55a 100644 --- a/man/ModelArray.gam.Rd +++ b/man/ModelArray.gam.Rd @@ -104,12 +104,13 @@ skips. Default: \code{"stop"}.} \item{write_results_name}{Optional character. If provided, results are incrementally written to -\code{results//results_matrix} in the HDF5 file +\code{results//results_matrix} in the output specified by \code{write_results_file}.} -\item{write_results_file}{Optional character. HDF5 file path for -incremental result writes. Required when \code{write_results_name} -is provided.} +\item{write_results_file}{Optional character. HDF5 file path or TileDB +store path for incremental result writes. Required when +\code{write_results_name} is provided. The backend is detected from the +path, so use a \code{.tdb} suffix to stream into a TileDB store.} \item{write_results_flush_every}{Positive integer. Number of elements per write block. Default 1000.} @@ -122,7 +123,7 @@ level for HDF5 writes. Default 4.} \item{return_output}{Logical. If \code{TRUE} (default), return the combined data.frame. If \code{FALSE}, return \code{invisible(NULL)}; -useful when writing large outputs directly to HDF5.} +useful when writing outputs directly to storage.} \item{...}{Additional arguments passed to \code{\link[mgcv]{gam}} (e.g. \code{method = "REML"}).} diff --git a/man/ModelArray.lm.Rd b/man/ModelArray.lm.Rd index af6dd6d..fd706ad 100644 --- a/man/ModelArray.lm.Rd +++ b/man/ModelArray.lm.Rd @@ -90,12 +90,13 @@ skips. Default: \code{"stop"}.} \item{write_results_name}{Optional character. If provided, results are incrementally written to -\code{results//results_matrix} in the HDF5 file +\code{results//results_matrix} in the output specified by \code{write_results_file}.} -\item{write_results_file}{Optional character. HDF5 file path for -incremental result writes. Required when \code{write_results_name} -is provided.} +\item{write_results_file}{Optional character. HDF5 file path or TileDB +store path for incremental result writes. Required when +\code{write_results_name} is provided. The backend is detected from the +path, so use a \code{.tdb} suffix to stream into a TileDB store.} \item{write_results_flush_every}{Positive integer. Number of elements per write block. Default 1000.} @@ -108,7 +109,7 @@ level for HDF5 writes. Default 4.} \item{return_output}{Logical. If \code{TRUE} (default), return the combined data.frame. If \code{FALSE}, return \code{invisible(NULL)}; -useful when writing large outputs directly to HDF5.} +useful when writing outputs directly to storage.} \item{...}{Additional arguments passed to \code{\link[stats]{lm}}.} } diff --git a/man/ModelArray.wrap.Rd b/man/ModelArray.wrap.Rd index 79e5103..bf63efa 100644 --- a/man/ModelArray.wrap.Rd +++ b/man/ModelArray.wrap.Rd @@ -76,11 +76,13 @@ skips. Default: \code{"stop"}.} \item{write_scalar_name}{Optional character. If provided, selected output columns are written into -\code{scalars//values} in the HDF5 file specified -by \code{write_scalar_file}.} +\code{scalars//values} in the output specified by +\code{write_scalar_file}.} -\item{write_scalar_file}{Optional character. HDF5 output file path. -Required when \code{write_scalar_name} is provided.} +\item{write_scalar_file}{Optional character. HDF5 file path or TileDB store +path. Required when \code{write_scalar_name} is provided. The backend is +detected from the path, so use a \code{.tdb} suffix to write a TileDB +store.} \item{write_scalar_columns}{Optional character or integer vector selecting which output columns to save as scalar values. If @@ -95,19 +97,21 @@ names saved as scalar \code{column_names}. If \code{NULL}, uses block for scalar writes. Default 1000.} \item{write_scalar_storage_mode}{Character. Storage mode for scalar -writes (e.g. \code{"double"}). Default \code{"double"}.} +writes (e.g. \code{"double"}). TileDB supports \code{"logical"}, +\code{"integer"}, and \code{"double"}. Default \code{"double"}.} \item{write_scalar_compression_level}{Integer 0--9. Gzip compression -level for scalar writes. Default 4.} +level for HDF5 and TileDB scalar writes. Default 4.} \item{write_results_name}{Optional character. If provided, results are incrementally written to -\code{results//results_matrix} in the HDF5 file +\code{results//results_matrix} in the output specified by \code{write_results_file}.} -\item{write_results_file}{Optional character. HDF5 file path for -incremental result writes. Required when \code{write_results_name} -is provided.} +\item{write_results_file}{Optional character. HDF5 file path or TileDB +store path for incremental result writes. Required when +\code{write_results_name} is provided. The backend is detected from the +path, so use a \code{.tdb} suffix to stream into a TileDB store.} \item{write_results_flush_every}{Positive integer. Number of elements per write block. Default 1000.} @@ -120,7 +124,7 @@ level for HDF5 writes. Default 4.} \item{return_output}{Logical. If \code{TRUE} (default), return the combined data.frame. If \code{FALSE}, return \code{invisible(NULL)}; -useful when writing large outputs directly to HDF5.} +useful when writing outputs directly to storage.} \item{...}{Additional arguments forwarded to \code{FUN}.} } diff --git a/man/ModelArraySummary.Rd b/man/ModelArraySummary.Rd new file mode 100644 index 0000000..f1d073d --- /dev/null +++ b/man/ModelArraySummary.Rd @@ -0,0 +1,80 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{ModelArraySummary} +\alias{ModelArraySummary} +\alias{h5summary} +\alias{print.ModelArraySummary} +\alias{print.h5summary} +\title{Summarize ModelArray storage without loading a full ModelArray} +\usage{ +ModelArraySummary(filepath, backend = c("auto", "hdf5", "tiledb")) + +h5summary(filepath, backend = c("auto", "hdf5", "tiledb")) + +\method{print}{ModelArraySummary}(x, ...) + +\method{print}{h5summary}(x, ...) +} +\arguments{ +\item{filepath}{Character. Path to an HDF5 (\code{.h5}) file or TileDB +(\code{.tdb}) store.} + +\item{backend}{Character. Storage backend: \code{"auto"} (default), +\code{"hdf5"}, or \code{"tiledb"}.} + +\item{x}{A \code{ModelArraySummary} or \code{h5summary} object.} + +\item{...}{Additional arguments (currently ignored).} +} +\value{ +\code{ModelArraySummary()} returns an object of class +\code{c("ModelArraySummary", "h5summary")}. \code{h5summary()} returns +an object of class \code{"h5summary"} for backward compatibility. Both +are lists with components: +\describe{ +\item{scalars}{A data.frame with columns \code{name}, +\code{nElements}, and \code{nInputFiles}.} +\item{analyses}{Character vector of analysis names found under +\code{/results/}.} +\item{filepath}{The input filepath.} +\item{backend}{The resolved storage backend.} +} + +Invisible \code{x}. Called for its side effect of printing a +human-readable summary to the console. +} +\description{ +Reads the storage structure and returns a summary of available scalars, +their dimensions, and any saved analyses. Useful for inspecting large files +without constructing a full \linkS4class{ModelArray} object. +} +\details{ +For HDF5, this function opens the file read-only via +\code{\link[rhdf5]{h5ls}}. For TileDB, it inspects arrays under +\code{/scalars/} and \code{/results/}. It does not load full data matrices +into memory. The returned +object has a \code{print} method that displays a formatted summary. + +\code{h5summary()} is a backward-compatible alias for +\code{ModelArraySummary()}. Despite the historical name, it supports both +HDF5 files and TileDB stores. +} +\examples{ +\dontrun{ +ModelArraySummary("path/to/data.h5") +ModelArraySummary("path/to/store.tdb") + +# Inspect before deciding which scalars to load +info <- ModelArraySummary("path/to/data.h5") +info$scalars$name +ma <- ModelArray("path/to/data.h5", scalar_types = info$scalars$name) + +# Historical alias, still supported +h5summary("path/to/data.h5") +} + +} +\seealso{ +\code{\link{ModelArray}} for loading the full object, +\linkS4class{ModelArray} for the class definition. +} diff --git a/man/elementMetadata.Rd b/man/elementMetadata.Rd index 7590f4e..96932d9 100644 --- a/man/elementMetadata.Rd +++ b/man/elementMetadata.Rd @@ -14,13 +14,15 @@ elementMetadata(x) } \value{ A matrix or data.frame of element metadata if found, or \code{NULL} -if no known metadata dataset exists in the HDF5 file. +if no known metadata dataset exists. } \description{ Reads element metadata (e.g., greyordinates for CIFTI data, or fixel/voxel -coordinate information) from the HDF5 file if present. The function searches -for known metadata dataset names (\code{"greyordinates"}, \code{"fixels"}, -\code{"voxels"}) at the top level of the HDF5 file. +coordinate information) from the HDF5 file or TileDB store if present. The +function searches for known metadata dataset names (\code{"greyordinates"}, +\code{"fixels"}, \code{"voxels"}) at the top level of the storage backend. +It uses the backend retained when the \code{ModelArray} was constructed, so +an explicit constructor choice is preserved. } \examples{ \dontrun{ diff --git a/man/h5summary.Rd b/man/h5summary.Rd deleted file mode 100644 index 1788988..0000000 --- a/man/h5summary.Rd +++ /dev/null @@ -1,59 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R -\name{h5summary} -\alias{h5summary} -\alias{print.h5summary} -\title{Summarize an HDF5 file without loading a full ModelArray} -\usage{ -h5summary(filepath) - -\method{print}{h5summary}(x, ...) -} -\arguments{ -\item{filepath}{Character. Path to an HDF5 (\code{.h5}) file.} - -\item{x}{An \code{h5summary} object as returned by -\code{\link{h5summary}}.} - -\item{...}{Additional arguments (currently ignored).} -} -\value{ -An object of class \code{"h5summary"}, which is a list with -components: -\describe{ -\item{scalars}{A data.frame with columns \code{name}, -\code{nElements}, and \code{nInputFiles}.} -\item{analyses}{Character vector of analysis names found under -\code{/results/}.} -\item{filepath}{The input filepath.} -} - -Invisible \code{x}. Called for its side effect of printing a -human-readable summary to the console. -} -\description{ -Reads the HDF5 file structure and returns a summary of available scalars, -their dimensions, and any saved analyses. Useful for inspecting large files -without constructing a full \linkS4class{ModelArray} object. -} -\details{ -This function opens the HDF5 file read-only via \code{\link[rhdf5]{h5ls}}, -inspects the group structure under \code{/scalars/} and \code{/results/}, -and closes the file. It does not load any data into memory. The returned -object has a \code{print} method that displays a formatted summary. -} -\examples{ -\dontrun{ -h5summary("path/to/data.h5") - -# Inspect before deciding which scalars to load -info <- h5summary("path/to/data.h5") -info$scalars$name -ma <- ModelArray("path/to/data.h5", scalar_types = info$scalars$name) -} - -} -\seealso{ -\code{\link{ModelArray}} for loading the full object, -\linkS4class{ModelArray} for the class definition. -} diff --git a/man/writeResults.Rd b/man/writeResults.Rd index 9e0653e..21fd749 100644 --- a/man/writeResults.Rd +++ b/man/writeResults.Rd @@ -2,19 +2,20 @@ % Please edit documentation in R/ModelArray_Constructor.R \name{writeResults} \alias{writeResults} -\title{Write outputs from element-wise statistical analysis to an HDF5 file} +\title{Write outputs from element-wise statistical analysis to storage} \usage{ writeResults( fn.output, df.output, analysis_name = "myAnalysis", - overwrite = TRUE + overwrite = TRUE, + backend = c("auto", "hdf5", "tiledb") ) } \arguments{ -\item{fn.output}{Character. The HDF5 (\code{.h5}) filename for the output. -The file must already exist; use an absolute path if you encounter -file-not-found errors.} +\item{fn.output}{Character. The HDF5 (\code{.h5}) filename or TileDB +(\code{.tdb}) store for the output. Use an absolute path if you +encounter file-not-found errors.} \item{df.output}{A data.frame of element-wise statistical results, as returned by \code{\link{ModelArray.lm}}, @@ -26,23 +27,29 @@ as the group name under \code{/results/} in the HDF5 file. Default is \code{"myAnalysis"}.} \item{overwrite}{Logical. If a group with the same \code{analysis_name} -already exists in the HDF5 file, whether to overwrite it (\code{TRUE}) +already exists in the storage backend, whether to overwrite it (\code{TRUE}) or skip with a warning (\code{FALSE}). Default is \code{TRUE}.} + +\item{backend}{Character. Storage backend: \code{"auto"} (default), +\code{"hdf5"}, or \code{"tiledb"}. Auto-detection resolves to TileDB for +a \code{.tdb} path, or for a directory that already contains a +\code{scalars/} or \code{results/} subdirectory; everything else is +treated as HDF5. To create a new TileDB store at a path without a +\code{.tdb} suffix, pass \code{"tiledb"} explicitly.} } \value{ Invisible \code{NULL}. Called for its side effect of writing -results to the HDF5 file. +results. } \description{ Creates a group named \code{analysis_name} under \code{/results/} in the -HDF5 file, then writes the statistical results data.frame (i.e. for one +storage backend, then writes the statistical results data.frame (i.e. for one analysis) into it as \code{results_matrix} along with column names. } \details{ The results are stored at \code{/results//results_matrix} with column names saved -as a separate dataset at -\code{/results//column_names}. +as a separate HDF5 dataset or TileDB metadata. If any column of \code{df.output} is not numeric or integer, it is coerced to numeric via \code{factor()} and the factor levels are saved @@ -74,7 +81,7 @@ writeResults( ) # Verify -h5summary("data.h5") +ModelArraySummary("data.h5") } } @@ -82,6 +89,6 @@ h5summary("data.h5") \code{\link{ModelArray.lm}}, \code{\link{ModelArray.gam}}, \code{\link{ModelArray.wrap}} which produce the \code{df.output}, \code{\link{results}} for reading results back from a -\linkS4class{ModelArray}, \code{\link{h5summary}} for inspecting what -has been written. +\linkS4class{ModelArray}, \code{\link{ModelArraySummary}} for inspecting +what has been written. } diff --git a/tests/testthat/helper-expected_values.R b/tests/testthat/helper-expected_values.R index 49af270..456f36f 100644 --- a/tests/testthat/helper-expected_values.R +++ b/tests/testthat/helper-expected_values.R @@ -76,7 +76,7 @@ calcu_stat_gam <- function(formula, data, i_element = idx.fixel.gam, ...) { onemodel.tidy.parametricTerms <- onemodel %>% broom::tidy(parametric = TRUE) # needs to use broom::glance instead of broom::glance.gam()! onemodel.glance <- onemodel %>% broom::glance() - onemodel.summary <- onemodel %>% summary.gam() + onemodel.summary <- onemodel %>% mgcv::summary.gam() # add additional model's stat to onemodel.glance(): onemodel.glance[["adj.r.squared"]] <- onemodel.summary$r.sq onemodel.glance[["dev.expl"]] <- onemodel.summary$dev.expl @@ -142,17 +142,17 @@ calcu_stat_gam <- function(formula, data, i_element = idx.fixel.gam, ...) { # all(union) = c(TRUE, FALSE, ...) temp <- union(temp_colnames, "term") # just an empty tibble (so below, all(dim(onemodel.tidy.smoothTerms)) = FALSE) - if (all(temp == "term")) onemodel.tidy.smoothTerms <- tibble() + if (all(temp == "term")) onemodel.tidy.smoothTerms <- tibble::tibble() temp_colnames <- onemodel.tidy.parametricTerms %>% colnames() temp <- union(temp_colnames, "term") # just an empty tibble - if (all(temp == "term")) onemodel.tidy.parametricTerms <- tibble() + if (all(temp == "term")) onemodel.tidy.parametricTerms <- tibble::tibble() temp_colnames <- onemodel.glance %>% colnames() temp <- union(temp_colnames, "term") # just an empty tibble - if (all(temp == "term")) onemodel.glance <- tibble() + if (all(temp == "term")) onemodel.glance <- tibble::tibble() ## flatten: @@ -256,7 +256,7 @@ helper_generate_expect_lm <- function(fn.phenotypes, fd.simu <- h5d[idx.fixel.lm, ] - h5closeAll() + rhdf5::h5closeAll() data <- phenotypes @@ -355,7 +355,7 @@ helper_generate_expect_gam <- function(fn.phenotypes, fd.simu <- h5d[idx.fixel.gam, ] - h5closeAll() + rhdf5::h5closeAll() data <- phenotypes data$FD <- fd.simu @@ -554,7 +554,7 @@ helper_gen_exp_accessors <- function(fn.h5) { scalar.value.full <- h5f$scalars$FD$values # enter the dataset # ^ has been realized as matrix, taking ~70MB of memory if using n50_fixels.h5 file - h5closeAll() + rhdf5::h5closeAll() scalar.value.full } diff --git a/tests/testthat/test-ModelArray_class.R b/tests/testthat/test-ModelArray_class.R index bfb77b8..74a91b1 100644 --- a/tests/testthat/test-ModelArray_class.R +++ b/tests/testthat/test-ModelArray_class.R @@ -11,7 +11,7 @@ test_that("ModelArray interface works as expected", { # ^ has been realized as matrix, taking ~70MB of memory if using n50_fixels.h5 file }, finally = { # regardless try is successful or not: - h5closeAll() + rhdf5::h5closeAll() } ) diff --git a/tests/testthat/test-ModelArray_gam.R b/tests/testthat/test-ModelArray_gam.R index be9f1f2..cb9f1db 100644 --- a/tests/testthat/test-ModelArray_gam.R +++ b/tests/testthat/test-ModelArray_gam.R @@ -99,10 +99,10 @@ test_that("test that ModelArray.gam() works as expected", { } # now, remove any fdr, etc corrections: - actual <- actual %>% select(-ends_with(p.adjust.methods)) + actual <- actual %>% dplyr::select(-ends_with(p.adjust.methods)) # also remove any effect size calculations: - actual <- actual %>% select(-ends_with(c( + actual <- actual %>% dplyr::select(-ends_with(c( "delta.adj.rsq", "partial.rsq" ))) @@ -120,7 +120,7 @@ test_that("test that ModelArray.gam() works as expected", { ## test if actual = expected values: expect_equal( actual, - expected %>% select(col.names) + expected %>% dplyr::select(col.names) ) } @@ -658,8 +658,8 @@ test_that("test that ModelArray.gam() works as expected", { mygam_fullOutputs %>% colnames() ) expect_equal( - mygam_changedrsq_oneSmoothTerm_sex %>% select(colnames_intersect), - mygam_fullOutputs %>% select(colnames_intersect) + mygam_changedrsq_oneSmoothTerm_sex %>% dplyr::select(colnames_intersect), + mygam_fullOutputs %>% dplyr::select(colnames_intersect) ) # compared to fdr: expect_equal( @@ -750,8 +750,8 @@ test_that("test that ModelArray.gam() works as expected", { mygam_twoSmoothTerm_withoutChangedRsq %>% colnames() ) expect_equal( - mygam_changedRsq_twoSmoothTerm %>% select(colnames_intersect), - mygam_twoSmoothTerm_withoutChangedRsq %>% select(colnames_intersect) + mygam_changedRsq_twoSmoothTerm %>% dplyr::select(colnames_intersect), + mygam_twoSmoothTerm_withoutChangedRsq %>% dplyr::select(colnames_intersect) ) # compared to fdr: expect_equal( diff --git a/tests/testthat/test-ModelArray_low_hanging_coverage.R b/tests/testthat/test-ModelArray_low_hanging_coverage.R index c439e1e..e888a8e 100644 --- a/tests/testthat/test-ModelArray_low_hanging_coverage.R +++ b/tests/testthat/test-ModelArray_low_hanging_coverage.R @@ -20,9 +20,9 @@ test_that("S4 accessors and helpers cover low-hanging branches", { # show() formatting branch shown <- paste(capture.output(show(modelarray)), collapse = "\n") - expect_match(shown, "elements x") - expect_match(shown, "input files") - expect_match(shown, "Analyses:") + expect_match(shown, "Source files:\\s+2") + expect_match(shown, "Scalars:\\s+FD \\(3 elements\\), FA \\(3 elements\\)") + expect_match(shown, "Analyses:\\s+my_analysis") # helper success + error branches phenotypes <- data.frame(source_file = src, age = c(20, 30)) @@ -78,6 +78,25 @@ test_that("ModelArray constructor falls back to dataset column names and transpo expect_identical(results(modelarray), list()) }) +test_that("elementMetadata honors the constructor's explicit backend", { + h5_path <- tempfile(fileext = ".tdb") + on.exit(unlink(h5_path), add = TRUE) + + metadata <- matrix(c(10, 20, 30, 40, 50, 60), nrow = 3, byrow = TRUE) + h5 <- hdf5r::H5File$new(h5_path, mode = "w") + scalars_grp <- h5$create_group("scalars") + fd_grp <- scalars_grp$create_group("FD") + fd_grp[["values"]] <- matrix(1:6, nrow = 3, ncol = 2) + fd_grp[["column_names"]] <- c("subA", "subB") + h5[["fixels"]] <- metadata + h5$close_all() + + modelarray <- ModelArray(h5_path, scalar_types = "FD", backend = "hdf5") + + expect_identical(modelarray@backend, "hdf5") + expect_equal(elementMetadata(modelarray), metadata) +}) + test_that("ModelArray constructor errors when scalar column names are unavailable", { h5_path <- tempfile(fileext = ".h5") on.exit(unlink(h5_path), add = TRUE) @@ -102,22 +121,12 @@ test_that("ModelArray constructor supports backward-compatible column-name attri rhdf5::h5createGroup(h5_path, "scalars") rhdf5::h5createGroup(h5_path, "scalars/FD") rhdf5::h5write(matrix(1:6, nrow = 3, ncol = 2), h5_path, "scalars/FD/values") - rhdf5::h5writeAttribute( - attr = c("attr_sub1", "attr_sub2"), - h5obj = h5_path, - name = "column_names", - h5loc = "scalars/FD/values" - ) + .write_hdf5_attribute(h5_path, "scalars/FD/values", "column_names", c("attr_sub1", "attr_sub2")) rhdf5::h5createGroup(h5_path, "results") rhdf5::h5createGroup(h5_path, "results/my_analysis") rhdf5::h5write(matrix(c(10, 20, 30, 40), nrow = 2, ncol = 2), h5_path, "results/my_analysis/results_matrix") - rhdf5::h5writeAttribute( - attr = c("beta", "p.value"), - h5obj = h5_path, - name = "colnames", - h5loc = "results/my_analysis/results_matrix" - ) + .write_hdf5_attribute(h5_path, "results/my_analysis/results_matrix", "colnames", c("beta", "p.value")) modelarray <- ModelArray( h5_path, @@ -141,23 +150,18 @@ test_that("ModelArray constructor prefers attributes over dataset column-name fa rhdf5::h5createGroup(h5_path, "scalars") rhdf5::h5createGroup(h5_path, "scalars/FD") rhdf5::h5write(matrix(1:6, nrow = 3, ncol = 2), h5_path, "scalars/FD/values") - rhdf5::h5writeAttribute( - attr = c("attr_sub1", "attr_sub2"), - h5obj = h5_path, - name = "column_names", - h5loc = "scalars/FD/values" - ) + .write_hdf5_attribute(h5_path, "scalars/FD/values", "column_names", c("attr_sub1", "attr_sub2")) # Conflicting fallback path should be ignored when attribute exists. rhdf5::h5write(c("fallback_sub1", "fallback_sub2"), h5_path, "scalars/FD/column_names") rhdf5::h5createGroup(h5_path, "results") rhdf5::h5createGroup(h5_path, "results/my_analysis") rhdf5::h5write(matrix(c(10, 20, 30, 40), nrow = 2, ncol = 2), h5_path, "results/my_analysis/results_matrix") - rhdf5::h5writeAttribute( - attr = c("attr_beta", "attr_p.value"), - h5obj = h5_path, - name = "colnames", - h5loc = "results/my_analysis/results_matrix" + .write_hdf5_attribute( + h5_path, + "results/my_analysis/results_matrix", + "colnames", + c("attr_beta", "attr_p.value") ) # Conflicting fallback path should be ignored when attribute exists. rhdf5::h5write(c("fallback_beta", "fallback_p"), h5_path, "results/my_analysis/column_names") diff --git a/tests/testthat/test-analyse_context.R b/tests/testthat/test-analyse_context.R index 4b1430b..c5804e9 100644 --- a/tests/testthat/test-analyse_context.R +++ b/tests/testthat/test-analyse_context.R @@ -85,7 +85,7 @@ test_that(".build_base_context attaches only requested scalars when scalar_subse test_that(".build_base_context detects collision between scalar names and phenotype columns", { fix <- make_test_modelarray(scalars = c("FD")) phen_collision <- fix$phen - phen_collision$FD <- 1:nrow(phen_collision) + phen_collision$FD <- seq_len(nrow(phen_collision)) expect_error( .build_base_context(fix$ma, phen_collision, scalar = "FD", @@ -311,6 +311,187 @@ test_that(".assemble_element_data intersection mask works across scalars", { expect_false(any(is.nan(elem$dat$FC))) }) +test_that(".with_scalar_row_cache serves cached rows for TileDB-style stores", { + store <- tempfile(fileext = ".tdb") + dir.create(store) + on.exit(unlink(store, recursive = TRUE), add = TRUE) + + old_block_size <- getOption("ModelArray.tiledb_read_block_size") + on.exit(options(ModelArray.tiledb_read_block_size = old_block_size), add = TRUE) + options(ModelArray.tiledb_read_block_size = 2L) + + src <- paste0("sub", 1:4) + values <- matrix(seq_len(12), nrow = 3, ncol = 4) + colnames(values) <- src + ma <- methods::new("ModelArray", + sources = list(FD = src), + scalars = list(FD = values), + results = list(), + path = store + ) + phen <- data.frame( + source_file = src, + age = seq(20, by = 10, length.out = length(src)) + ) + + ctx <- .build_lm_context(FD ~ age, ma, phen, scalar = "FD") + cached_ctx <- .with_scalar_row_cache(ctx, as.integer(c(3, 1))) + + expect_named(cached_ctx$scalar_row_cache, "FD") + expect_equal(cached_ctx$scalar_row_cache[["FD"]], values[c(3, 1), , drop = FALSE]) + + cached_ctx$modelarray@scalars[["FD"]][, ] <- 999 + elem <- .assemble_element_data(3L, cached_ctx, num.subj.lthr = 0) + + expect_true(elem$sufficient) + expect_equal(as.numeric(elem$dat$FD), as.numeric(values[3, ])) +}) + +test_that("TileDB row caching honors the backend retained by ModelArray", { + path_with_tiledb_suffix <- tempfile(fileext = ".tdb") + file.create(path_with_tiledb_suffix) + on.exit(unlink(path_with_tiledb_suffix), add = TRUE) + + src <- c("sub1", "sub2") + values <- matrix(seq_len(4), nrow = 2, ncol = 2) + colnames(values) <- src + ma <- methods::new( + "ModelArray", + sources = list(FD = src), + scalars = list(FD = values), + results = list(), + path = path_with_tiledb_suffix, + backend = "hdf5" + ) + phen <- data.frame(source_file = src, age = c(20, 30)) + ctx <- .build_lm_context(FD ~ age, ma, phen, scalar = "FD") + + expect_identical(.modelarray_backend(ma, path_with_tiledb_suffix), "hdf5") + expect_false(.has_tiledb_attached_scalar(ctx)) + expect_null(.scalar_read_block_size(ctx)) +}) + +test_that("TileDB read block sizing is memory bounded without a hard row cap", { + store <- tempfile(fileext = ".tdb") + dir.create(store) + on.exit(unlink(store, recursive = TRUE), add = TRUE) + + old_block_size <- getOption("ModelArray.tiledb_read_block_size") + old_block_mb <- getOption("ModelArray.tiledb_read_block_mb") + on.exit(options( + ModelArray.tiledb_read_block_size = old_block_size, + ModelArray.tiledb_read_block_mb = old_block_mb + ), add = TRUE) + options( + ModelArray.tiledb_read_block_size = NULL, + ModelArray.tiledb_read_block_mb = 1 + ) + + src <- paste0("sub", 1:2) + values <- matrix(seq_len(4), nrow = 2, ncol = 2) + colnames(values) <- src + ma <- methods::new("ModelArray", + sources = list(FD = src), + scalars = list(FD = values), + results = list(), + path = store + ) + phen <- data.frame(source_file = src, age = c(20, 30)) + ctx <- .build_lm_context(FD ~ age, ma, phen, scalar = "FD") + + expect_equal(.scalar_read_block_size(ctx), 65536L) + + options(ModelArray.tiledb_read_block_mb = -1) + expect_error( + .scalar_read_block_size(ctx), + "ModelArray.tiledb_read_block_mb must be a positive single number" + ) +}) + +test_that(".report_iteration_progress labels chunked progress bars", { + out <- capture.output( + .report_iteration_progress( + iteration = 1L, + total_iterations = 3L, + chunk_start = 1L, + chunk_end = 10L, + total_elements = 25L, + pbar = TRUE + ), + type = "message" + ) + expect_match(paste(out, collapse = "\n"), "iteration 1/3") + expect_match(paste(out, collapse = "\n"), "elements 1-10 of 25") + + quiet <- capture.output( + .report_iteration_progress(1L, 1L, 1L, 10L, 10L, TRUE), + type = "message" + ) + expect_length(quiet, 0L) +}) + +test_that("ModelArray.lm returns same output with TileDB-style row cache enabled", { + store <- tempfile(fileext = ".tdb") + dir.create(store) + on.exit(unlink(store, recursive = TRUE), add = TRUE) + + old_block_size <- getOption("ModelArray.tiledb_read_block_size") + on.exit(options(ModelArray.tiledb_read_block_size = old_block_size), add = TRUE) + + src <- paste0("sub", 1:8) + values <- matrix(c( + 1, 3, 2, 7, 5, 9, 8, 6, + 2, 1, 4, 3, 6, 5, 7, 8, + 3, 2, 5, 4, 7, 6, 8, 9, + 4, 6, 5, 8, 9, 7, 10, 11 + ), nrow = 4, byrow = TRUE) + colnames(values) <- src + ma <- methods::new("ModelArray", + sources = list(FD = src), + scalars = list(FD = values), + results = list(), + path = store + ) + phen <- data.frame(source_file = src, age = seq(10, 80, by = 10)) + element_subset <- as.integer(c(4, 1, 3, 2)) + + options(ModelArray.tiledb_read_block_size = 2L) + cached <- ModelArray.lm( + FD ~ age, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + var.terms = c("estimate"), + var.model = c("adj.r.squared"), + correct.p.value.terms = c("none"), + correct.p.value.model = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE + ) + + options(ModelArray.tiledb_read_block_size = 0L) + uncached <- ModelArray.lm( + FD ~ age, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + var.terms = c("estimate"), + var.model = c("adj.r.squared"), + correct.p.value.terms = c("none"), + correct.p.value.model = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE + ) + + expect_equal(cached, uncached) +}) + # ================================================================== # SECTION 3: Per-element function ctx vs. legacy equivalence diff --git a/tests/testthat/test-convenience_accessors.R b/tests/testthat/test-convenience_accessors.R index 2b75f61..7871710 100644 --- a/tests/testthat/test-convenience_accessors.R +++ b/tests/testthat/test-convenience_accessors.R @@ -37,10 +37,11 @@ test_that("elementMetadata returns fixel metadata", { expect_equal(nrow(em), 182581L) }) -test_that("h5summary returns correct structure", { +test_that("ModelArraySummary returns correct structure", { h5_path <- system.file("extdata", "n50_fixels.h5", package = "ModelArray") - s <- h5summary(h5_path) + s <- ModelArraySummary(h5_path) + expect_s3_class(s, "ModelArraySummary") expect_s3_class(s, "h5summary") expect_equal(s$scalars$name, "FD") expect_equal(s$scalars$nElements, 182581L) @@ -53,6 +54,24 @@ test_that("h5summary returns correct structure", { expect_true(any(grepl("182581", output))) }) +test_that("h5summary remains a backward-compatible alias", { + h5_path <- system.file("extdata", "n50_fixels.h5", package = "ModelArray") + + current <- ModelArraySummary(h5_path) + legacy <- h5summary(h5_path) + + expect_s3_class(legacy, "h5summary") + expect_false(inherits(legacy, "ModelArraySummary")) + expect_equal(legacy$scalars, current$scalars) + expect_equal(legacy$analyses, current$analyses) + expect_equal(legacy$filepath, current$filepath) + expect_equal(legacy$backend, current$backend) +}) + test_that("h5summary errors on missing file", { expect_error(h5summary("/nonexistent/file.h5"), "not found") }) + +test_that("ModelArraySummary errors on missing file", { + expect_error(ModelArraySummary("/nonexistent/file.h5"), "not found") +}) diff --git a/tests/testthat/test-stream_results.R b/tests/testthat/test-stream_results.R index 92e38f0..6e9dfe3 100644 --- a/tests/testthat/test-stream_results.R +++ b/tests/testthat/test-stream_results.R @@ -122,3 +122,21 @@ test_that("ModelArray.wrap can stream results to HDF5 without returning output", mat_h5 <- rhdf5::h5read(h5_out, "results/wrap_stream/results_matrix") expect_equal(dim(mat_h5), c(2, 3)) }) + +test_that(".init_results_stream_writer validates compression level for both backends", { + # Validation is hoisted above the backend branch, so an out-of-range level is + # rejected for HDF5 too, not only for TileDB. + for (path in c(tempfile(fileext = ".h5"), tempfile(fileext = ".tdb"))) { + expect_error( + .init_results_stream_writer( + write_results_name = "an_analysis", + write_results_file = path, + n_rows = 2L, + column_names = c("element_id", "estimate"), + compression_level = 10L + ), + "write_results_compression_level must be a single integer between 0 and 9" + ) + expect_false(file.exists(path)) + } +}) diff --git a/tests/testthat/test-tiledb_backend.R b/tests/testthat/test-tiledb_backend.R new file mode 100644 index 0000000..7b5b6d5 --- /dev/null +++ b/tests/testthat/test-tiledb_backend.R @@ -0,0 +1,971 @@ +skip_if_no_tiledb <- function() { + testthat::skip_if_not_installed("tiledb") + testthat::skip_if_not_installed("TileDBArray") + testthat::skip_if_not_installed("jsonlite") +} + +write_tiledb_group_metadata <- function(uri, key, value) { + dir.create(uri, recursive = TRUE, showWarnings = FALSE) + if (!identical(toupper(tiledb::tiledb_object_type(uri)), "GROUP")) { + tiledb::tiledb_group_create(uri) + } + grp <- tiledb::tiledb_group(uri, type = "WRITE") + on.exit(try(tiledb::tiledb_group_close(grp), silent = TRUE), add = TRUE) + tiledb::tiledb_group_put_metadata( + grp, + key, + as.character(jsonlite::toJSON(value, auto_unbox = TRUE)) + ) +} + +write_tiledb_scalar <- function(store, scalar, values, sources, metadata_location = c("array", "group")) { + metadata_location <- match.arg(metadata_location) + dir.create(file.path(store, "scalars", scalar), recursive = TRUE, showWarnings = FALSE) + uri <- file.path(store, "scalars", scalar, "values") + TileDBArray::writeTileDBArray(values, path = uri, attr = "values") + + dir.create(file.path(store, "scalars", "__group"), showWarnings = FALSE) + dir.create(file.path(store, "scalars", "__meta"), showWarnings = FALSE) + + if (identical(metadata_location, "array")) { + arr <- tiledb::tiledb_array(uri, query_type = "WRITE", keep_open = TRUE) + on.exit(try(tiledb::tiledb_array_close(arr), silent = TRUE), add = TRUE) + tiledb::tiledb_put_metadata( + arr, + "column_names", + as.character(jsonlite::toJSON(sources, auto_unbox = TRUE)) + ) + } else { + write_tiledb_group_metadata(file.path(store, "scalars", scalar), "column_names", sources) + } + invisible(uri) +} + +write_tiledb_result_with_group_metadata <- function(store, analysis_name, matrix_values, column_names, luts = list()) { + analysis_uri <- file.path(store, "results", analysis_name) + matrix_uri <- file.path(analysis_uri, "results_matrix") + dir.create(analysis_uri, recursive = TRUE, showWarnings = FALSE) + TileDBArray::writeTileDBArray(matrix_values, path = matrix_uri, attr = "values") + write_tiledb_group_metadata(analysis_uri, "column_names", column_names) + for (lut_name in names(luts)) { + write_tiledb_group_metadata(analysis_uri, lut_name, luts[[lut_name]]) + } + invisible(matrix_uri) +} + +test_that("ModelArray reads TileDB scalar arrays written subjects by elements", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + values <- matrix(seq_len(8), nrow = 2, byrow = TRUE) + sources <- c("sub-1.nii.gz", "sub-2.nii.gz") + write_tiledb_scalar(store, "FD", values, sources) + + ma <- ModelArray(store, scalar_types = "FD") + + expect_equal(sources(ma)[["FD"]], sources) + loaded <- as.matrix(scalars(ma)[["FD"]]) + dimnames(loaded) <- NULL + expect_equal(loaded, t(values)) +}) + +test_that("ModelArraySummary summarizes TileDB scalar arrays", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + write_tiledb_scalar( + store, + "FA", + matrix(seq_len(12), nrow = 3), + paste0("sub-", 1:3) + ) + + info <- ModelArraySummary(store, backend = "tiledb") + + expect_s3_class(info, "ModelArraySummary") + expect_s3_class(info, "h5summary") + expect_equal(info$backend, "tiledb") + expect_equal(info$scalars$name, "FA") + expect_equal(info$scalars$nInputFiles, 3L) + expect_equal(info$scalars$nElements, 4L) + + legacy <- h5summary(store, backend = "tiledb") + expect_s3_class(legacy, "h5summary") + expect_false(inherits(legacy, "ModelArraySummary")) + expect_equal(legacy$scalars, info$scalars) + expect_equal(legacy$backend, "tiledb") +}) + +test_that("ModelArraySummary agrees for fake HDF5 and TileDB stores", { + skip_if_no_tiledb() + + h5_store <- tempfile(fileext = ".h5") + tiledb_store <- tempfile(fileext = ".tdb") + on.exit(unlink(h5_store), add = TRUE) + on.exit(unlink(tiledb_store, recursive = TRUE), add = TRUE) + + values <- matrix(seq_len(6), nrow = 3, ncol = 2) + source_names <- c("sub-1", "sub-2") + for (store in c(h5_store, tiledb_store)) { + writer <- .init_scalar_stream_writer( + write_scalar_name = "FA", + write_scalar_file = store, + n_rows = nrow(values), + column_names = source_names, + flush_every = 2L + ) + writer <- .scalar_stream_write_block(writer, values[1:2, , drop = FALSE]) + writer <- .scalar_stream_write_block(writer, values[3, , drop = FALSE]) + .finalize_scalar_stream_writer(writer) + + writeResults( + store, + data.frame(element_id = 0:2, estimate = c(0.1, 0.2, 0.3)), + analysis_name = "fake_analysis" + ) + } + + h5_info <- ModelArraySummary(h5_store, backend = "hdf5") + tiledb_info <- ModelArraySummary(tiledb_store, backend = "tiledb") + + expect_equal( + h5_info[c("scalars", "analyses")], + tiledb_info[c("scalars", "analyses")] + ) + expect_equal(h5_info$scalars$nElements, 3L) + expect_equal(h5_info$scalars$nInputFiles, 2L) +}) + +test_that("writeResults can round-trip TileDB results metadata and LUTs", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + write_tiledb_scalar( + store, + "FD", + matrix(c(1, 2, 3, 4), nrow = 2, byrow = TRUE), + c("s1", "s2") + ) + + df <- data.frame( + element_id = 0:1, + beta = c(0.25, 0.5), + gamma = c(1.25, 1.5), + label = c("low", "high") + ) + writeResults(store, df, analysis_name = "lm_tile", backend = "tiledb") + + ma <- ModelArray( + store, + scalar_types = "FD", + analysis_names = "lm_tile", + backend = "tiledb" + ) + result_matrix <- results(ma, "lm_tile")$results_matrix + expect_s4_class(result_matrix, "DelayedMatrix") + expect_false(is.matrix(result_matrix)) + loaded <- as.matrix(result_matrix) + + expect_equal(colnames(loaded), colnames(df)) + expect_equal(as.numeric(loaded[, "element_id"]), df$element_id) + expect_equal(as.numeric(loaded[, "beta"]), df$beta) + expect_equal(as.numeric(loaded[, "gamma"]), df$gamma) + expect_equal(loaded[, "label"], df$label) +}) + +test_that("TileDB writers reject unsafe names without deleting stored data", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + dir.create(file.path(store, "scalars"), recursive = TRUE) + dir.create(file.path(store, "results"), recursive = TRUE) + scalar_marker <- file.path(store, "scalars", "keep.txt") + result_marker <- file.path(store, "results", "keep.txt") + writeLines("keep", scalar_marker) + writeLines("keep", result_marker) + result_data <- data.frame(element_id = 0, estimate = 0.25) + + unsafe_names <- c("", ".", "..", "../scalars", "nested/name", "nested\\name") + for (unsafe_name in unsafe_names) { + expect_error( + writeResults( + store, + result_data, + analysis_name = unsafe_name, + backend = "tiledb" + ), + "analysis_name" + ) + expect_true(file.exists(scalar_marker)) + expect_true(file.exists(result_marker)) + } + + expect_error( + .init_results_stream_writer( + write_results_name = "../scalars", + write_results_file = store, + n_rows = 1L, + column_names = names(result_data) + ), + "write_results_name" + ) + expect_true(file.exists(scalar_marker)) + expect_true(file.exists(result_marker)) + + expect_error( + .init_scalar_stream_writer( + write_scalar_name = "../results", + write_scalar_file = store, + n_rows = 1L, + column_names = "estimate" + ), + "write_scalar_name" + ) + expect_true(file.exists(scalar_marker)) + expect_true(file.exists(result_marker)) +}) + +test_that("ModelArray reads TileDB scalar and result metadata from groups", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + write_tiledb_scalar( + store, + "FD", + matrix(c(1, 2, 3, 4), nrow = 2, byrow = TRUE), + c("s1", "s2"), + metadata_location = "group" + ) + write_tiledb_result_with_group_metadata( + store, + "group_meta", + matrix(c(0, 1, 1, 2), nrow = 2, byrow = TRUE), + c("element_id", "label"), + luts = list(lut_forcol2 = c("low", "high")) + ) + + ma <- ModelArray( + store, + scalar_types = "FD", + analysis_names = "group_meta", + backend = "tiledb" + ) + loaded <- as.matrix(results(ma, "group_meta")$results_matrix) + + expect_equal(sources(ma)[["FD"]], c("s1", "s2")) + expect_equal(colnames(loaded), c("element_id", "label")) + expect_equal(loaded[, "label"], c("low", "high")) +}) + +test_that("analysis functions can consume TileDB-backed scalar arrays", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + src <- paste0("s", 1:8) + write_tiledb_scalar( + store, + "FD", + matrix(c( + 1, 3, 2, 7, 5, 9, 8, 6, + 2, 1, 4, 3, 6, 5, 7, 8 + ), nrow = 2, byrow = TRUE), + src + ) + ma <- ModelArray(store, scalar_types = "FD", backend = "tiledb") + phen <- data.frame(source_file = src, age = seq(10, 80, by = 10)) + + lm_out <- ModelArray.lm( + FD ~ age, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = as.integer(c(1, 2)), + var.terms = c("estimate"), + var.model = c("adj.r.squared"), + correct.p.value.terms = c("none"), + correct.p.value.model = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE + ) + gam_out <- ModelArray.gam( + FD ~ age, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = as.integer(c(1, 2)), + var.smoothTerms = c(), + var.parametricTerms = c("estimate"), + var.model = c("dev.expl"), + correct.p.value.smoothTerms = c("none"), + correct.p.value.parametricTerms = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE + ) + wrap_out <- ModelArray.wrap( + FUN = function(data) data.frame(mean_fd = mean(data$FD)), + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = as.integer(c(1, 2)), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE + ) + + expect_equal(nrow(lm_out), 2) + expect_equal(nrow(gam_out), 2) + expect_equal(wrap_out$mean_fd, rowMeans(as.matrix(scalars(ma, "FD")))) +}) + +test_that("TileDB result streaming writes each block without retaining it", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + writer <- .init_results_stream_writer( + write_results_name = "direct_stream", + write_results_file = store, + n_rows = 3L, + column_names = c("element_id", "estimate"), + flush_every = 1L + ) + + expect_false("blocks" %in% names(writer)) + writer <- .results_stream_write_block( + writer, + data.frame(element_id = 0, estimate = 0.25) + ) + expect_false("blocks" %in% names(writer)) + + fragment_info <- tiledb::tiledb_fragment_info(writer$values_uri) + expect_equal(tiledb::tiledb_fragment_info_get_num(fragment_info), 1L) + + writer <- .results_stream_write_block( + writer, + data.frame( + element_id = c(1, 2), + estimate = c(0.5, NaN) + ) + ) + .finalize_results_stream_writer(writer) + + fragment_info <- tiledb::tiledb_fragment_info(writer$values_uri) + expect_equal(tiledb::tiledb_fragment_info_get_num(fragment_info), 2L) +}) + +test_that("TileDB streamed outputs exactly match regular lm, gam, and wrap outputs", { + skip_if_no_tiledb() + + src <- paste0("s", 1:8) + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + write_tiledb_scalar( + store, + "FD", + matrix(c( + 1, 3, 2, 7, 5, 9, 8, 6, + 2, 1, 4, 3, 6, 5, 7, 8, + rep(NA_real_, 8) + ), nrow = 3, byrow = TRUE), + src + ) + ma <- ModelArray(store, scalar_types = "FD", backend = "tiledb") + phen <- data.frame(source_file = src, age = seq(10, 80, by = 10)) + element_subset <- as.integer(c(3, 1, 2)) + + expected_lm <- ModelArray.lm( + FD ~ age, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + var.terms = c("estimate", "std.error", "statistic", "p.value"), + var.model = c("adj.r.squared", "AIC"), + correct.p.value.terms = c("none"), + correct.p.value.model = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE + ) + + out_lm <- ModelArray.lm( + FD ~ age, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + var.terms = c("estimate", "std.error", "statistic", "p.value"), + var.model = c("adj.r.squared", "AIC"), + correct.p.value.terms = c("none"), + correct.p.value.model = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE, + write_results_name = "lm_stream", + write_results_file = store, + write_results_flush_every = 1L, + return_output = FALSE + ) + + expected_gam <- ModelArray.gam( + FD ~ age, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + var.smoothTerms = c(), + var.parametricTerms = c("estimate", "std.error", "statistic", "p.value"), + var.model = c("adj.r.squared", "dev.expl", "AIC"), + correct.p.value.smoothTerms = c("none"), + correct.p.value.parametricTerms = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE + ) + out_gam <- ModelArray.gam( + FD ~ age, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + var.smoothTerms = c(), + var.parametricTerms = c("estimate", "std.error", "statistic", "p.value"), + var.model = c("adj.r.squared", "dev.expl", "AIC"), + correct.p.value.smoothTerms = c("none"), + correct.p.value.parametricTerms = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE, + write_results_name = "gam_stream", + write_results_file = store, + write_results_flush_every = 1L, + return_output = FALSE + ) + + wrap_fun <- function(data) { + data.frame( + mean_fd = mean(data$FD), + sd_fd = stats::sd(data$FD) + ) + } + expected_wrap <- ModelArray.wrap( + FUN = wrap_fun, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE + ) + out_wrap <- ModelArray.wrap( + FUN = wrap_fun, + data = ma, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE, + write_results_name = "wrap_stream", + write_results_file = store, + write_results_flush_every = 1L, + return_output = FALSE + ) + + ma_with_results <- ModelArray( + store, + scalar_types = "FD", + analysis_names = c("lm_stream", "gam_stream", "wrap_stream"), + backend = "tiledb" + ) + expected_outputs <- list( + lm_stream = expected_lm, + gam_stream = expected_gam, + wrap_stream = expected_wrap + ) + for (analysis_name in names(expected_outputs)) { + saved <- as.matrix(results(ma_with_results, analysis_name)$results_matrix) + expected <- as.matrix(expected_outputs[[analysis_name]]) + expect_identical(colnames(saved), colnames(expected)) + expect_identical(unname(saved), unname(expected)) + + matrix_uri <- file.path( + store, + "results", + analysis_name, + "results_matrix" + ) + fragment_info <- tiledb::tiledb_fragment_info(matrix_uri) + expect_equal( + tiledb::tiledb_fragment_info_get_num(fragment_info), + length(element_subset) + ) + } + expect_null(out_lm) + expect_null(out_gam) + expect_null(out_wrap) +}) + +test_that("HDF5 and TileDB outputs are identical for regular and streamed writes", { + skip_if_no_tiledb() + + src <- paste0("s", 1:8) + values <- matrix(c( + 1, 3, 2, 7, 5, 9, 8, 6, + 2, 1, 4, 3, 6, 5, 7, 8, + rep(NA_real_, 8) + ), nrow = 3, byrow = TRUE) + modelarray <- methods::new( + "ModelArray", + sources = list(FD = src), + scalars = list(FD = values), + results = list(), + path = tempfile(fileext = ".h5") + ) + phen <- data.frame(source_file = src, age = seq(10, 80, by = 10)) + element_subset <- as.integer(c(3, 1, 2)) + + regular <- ModelArray.lm( + FD ~ age, + data = modelarray, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + var.terms = c("estimate", "std.error", "statistic", "p.value"), + var.model = c("adj.r.squared", "AIC"), + correct.p.value.terms = c("none"), + correct.p.value.model = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE + ) + + h5_regular_file <- tempfile(fileext = ".h5") + tiledb_regular_file <- tempfile(fileext = ".tdb") + h5_stream_file <- tempfile(fileext = ".h5") + tiledb_stream_file <- tempfile(fileext = ".tdb") + output_paths <- c( + h5_regular_file, + tiledb_regular_file, + h5_stream_file, + tiledb_stream_file + ) + on.exit(unlink(output_paths, recursive = TRUE), add = TRUE) + + writeResults( + h5_regular_file, + regular, + analysis_name = "regular", + backend = "hdf5" + ) + writeResults( + tiledb_regular_file, + regular, + analysis_name = "regular", + backend = "tiledb" + ) + + stream_lm <- function(output_file) { + ModelArray.lm( + FD ~ age, + data = modelarray, + phenotypes = phen, + scalar = "FD", + element.subset = element_subset, + var.terms = c("estimate", "std.error", "statistic", "p.value"), + var.model = c("adj.r.squared", "AIC"), + correct.p.value.terms = c("none"), + correct.p.value.model = c("none"), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE, + write_results_name = "stream", + write_results_file = output_file, + write_results_flush_every = 1L, + return_output = FALSE + ) + } + expect_null(stream_lm(h5_stream_file)) + expect_null(stream_lm(tiledb_stream_file)) + + read_results_matrix <- function(output_file, analysis_name, backend) { + saved <- as.matrix(DelayedArray::DelayedArray(ModelArraySeed( + output_file, + name = file.path("results", analysis_name, "results_matrix"), + backend = backend + ))) + column_names <- .read_result_column_names( + output_file, + analysis_name, + backend + ) + # Assert the stored orientation rather than normalizing it — silently + # transposing here would mask a backend writing elements x stats the + # wrong way round, which is the whole point of this comparison. + expect_identical( + ncol(saved), length(column_names), + info = paste("unexpected orientation for", backend, analysis_name) + ) + colnames(saved) <- column_names + saved + } + + h5_regular <- read_results_matrix(h5_regular_file, "regular", "hdf5") + tiledb_regular <- read_results_matrix( + tiledb_regular_file, + "regular", + "tiledb" + ) + h5_stream <- read_results_matrix(h5_stream_file, "stream", "hdf5") + tiledb_stream <- read_results_matrix( + tiledb_stream_file, + "stream", + "tiledb" + ) + + expect_identical(h5_regular, tiledb_regular) + expect_identical(h5_stream, tiledb_stream) + expect_identical(h5_regular, h5_stream) + expect_identical(tiledb_regular, tiledb_stream) +}) + +test_that("write_scalar streaming saves compressed scalars to TileDB", { + skip_if_no_tiledb() + + src <- paste0("s", 1:8) + modelarray <- methods::new( + "ModelArray", + sources = list(FD = src), + scalars = list(FD = matrix(c( + 1, 3, 2, 7, 5, 9, 8, 6, + 2, 1, 4, 3, 6, 5, 7, 8 + ), nrow = 2, byrow = TRUE)), + results = list(), + path = tempfile(fileext = ".h5") + ) + phen <- data.frame(source_file = src, age = seq(10, 80, by = 10)) + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + + out <- ModelArray.wrap( + FUN = function(data) data.frame(mean_fd = mean(data$FD)), + data = modelarray, + phenotypes = phen, + scalar = "FD", + element.subset = as.integer(c(1, 2)), + num.subj.lthr.abs = 0, + n_cores = 1, + pbar = FALSE, + verbose = FALSE, + write_scalar_name = "FD_harmonized", + write_scalar_file = store, + write_scalar_flush_every = 1L, + write_scalar_compression_level = 6L, + return_output = FALSE + ) + + saved_modelarray <- ModelArray( + store, + scalar_types = "FD_harmonized", + backend = "tiledb" + ) + saved <- as.matrix(scalars(saved_modelarray, "FD_harmonized")) + values_uri <- file.path(store, "scalars", "FD_harmonized", "values") + array <- tiledb::tiledb_array(values_uri, query_type = "READ", keep_open = TRUE) + on.exit(try(tiledb::tiledb_array_close(array), silent = TRUE), add = TRUE) + value_filters <- tiledb::filter_list( + tiledb::attrs(tiledb::schema(array))[["values"]] + ) + gzip_filter <- value_filters[0] + fragment_info <- tiledb::tiledb_fragment_info(values_uri) + + expect_null(out) + expect_equal(colnames(saved), "mean_fd") + expect_equal(as.vector(saved), rowMeans(modelarray@scalars$FD)) + expect_equal(tiledb::tiledb_fragment_info_get_num(fragment_info), 2L) + expect_equal(tiledb::nfilters(value_filters), 1L) + expect_equal(tiledb::tiledb_filter_type(gzip_filter), "GZIP") + expect_equal( + tiledb::tiledb_filter_get_option(gzip_filter, "COMPRESSION_LEVEL"), + 6L + ) +}) + +test_that("elementMetadata reads TileDB top-level metadata arrays", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + dir.create(store, recursive = TRUE) + metadata <- matrix(c(10, 20, 30, 40, 50, 60), nrow = 3, byrow = TRUE) + TileDBArray::writeTileDBArray( + metadata, + path = file.path(store, "fixels"), + attr = "values" + ) + ma <- methods::new( + "ModelArray", + sources = list(FD = c("s1", "s2")), + scalars = list(FD = matrix(1, nrow = 3, ncol = 2)), + results = list(), + path = store + ) + + expect_equal(elementMetadata(ma), metadata) +}) + +test_that("elementMetadata returns NULL when no TileDB metadata arrays exist", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + dir.create(store, recursive = TRUE) + ma <- methods::new( + "ModelArray", + sources = list(FD = c("s1", "s2")), + scalars = list(FD = matrix(1, nrow = 2, ncol = 2)), + results = list(), + path = store + ) + + expect_null(elementMetadata(ma)) +}) + +test_that("HDF5 and TileDB writeResults agree on classed-numeric columns", { + skip_if_no_tiledb() + + h5_file <- tempfile(fileext = ".h5") + store <- tempfile(fileext = ".tdb") + on.exit(unlink(h5_file), add = TRUE) + on.exit(unlink(store, recursive = TRUE), add = TRUE) + rhdf5::h5createFile(h5_file) + + # `scaled` is numeric underneath but carries an extra class attribute, so a + # deparsed-class check mistakes it for a factor candidate and replaces the + # values with LUT indices. Both backends must store the numbers as-is. + df <- data.frame(element_id = 0:2, estimate = c(0.5, 1.5, 2.5)) + df$scaled <- structure(c(10.5, 20.5, 30.5), class = c("myscale", "numeric")) + df$label <- c("a", "b", "a") + + writeResults(h5_file, df, analysis_name = "cmp", overwrite = TRUE) + writeResults(store, df, analysis_name = "cmp", overwrite = TRUE) + + read_back <- function(path, backend) { + ma <- ModelArray(path, + scalar_types = character(0), + analysis_names = "cmp", backend = backend + ) + as.matrix(results(ma)[["cmp"]]$results_matrix) + } + from_h5 <- read_back(h5_file, "hdf5") + from_tiledb <- read_back(store, "tiledb") + + expect_identical(from_h5, from_tiledb) + # numeric column survives intact rather than becoming 1/2/3 LUT indices + expect_equal(as.numeric(from_h5[, 3]), c(10.5, 20.5, 30.5)) + # the genuinely character column still round-trips through its LUT + expect_identical(as.character(from_h5[, 4]), df$label) + + # ...and the LUTs themselves match across backends + luts <- function(path, backend) { + vapply(seq_len(ncol(df)), function(i) { + lut <- .read_result_lut(path, "cmp", paste0("lut_forcol", i), backend) + if (is.null(lut)) NA_character_ else paste(lut, collapse = "|") + }, character(1)) + } + expect_identical(luts(h5_file, "hdf5"), luts(store, "tiledb")) +}) + +test_that("writeResults honors overwrite = FALSE on a TileDB store", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + + first <- data.frame(element_id = 0:1, estimate = c(1, 2)) + second <- data.frame(element_id = 0:1, estimate = c(99, 99)) + writeResults(store, first, analysis_name = "keep", overwrite = TRUE) + + expect_warning( + writeResults(store, second, analysis_name = "keep", overwrite = FALSE), + "exists but not to overwrite" + ) + kept <- ModelArray(store, scalar_types = character(0), analysis_names = "keep") + expect_equal( + as.numeric(as.matrix(results(kept)[["keep"]]$results_matrix)[, 2]), + first$estimate + ) + + writeResults(store, second, analysis_name = "keep", overwrite = TRUE) + replaced <- ModelArray(store, scalar_types = character(0), analysis_names = "keep") + expect_equal( + as.numeric(as.matrix(results(replaced)[["keep"]]$results_matrix)[, 2]), + second$estimate + ) +}) + +test_that("multi-scalar TileDB stores drive scalar predictors in a model", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + + src <- paste0("sub-", 1:6) + set.seed(7) + fd <- matrix(stats::rnorm(18), nrow = 3, ncol = 6) + fc <- matrix(stats::rnorm(18), nrow = 3, ncol = 6) + write_tiledb_scalar(store, "FD", fd, src) + write_tiledb_scalar(store, "FC", fc, src) + + ma <- ModelArray(store, scalar_types = c("FD", "FC")) + expect_identical(sort(names(scalars(ma))), c("FC", "FD")) + + phen <- data.frame(source_file = src, age = seq(10, 60, by = 10)) + out <- ModelArray.lm( + FD ~ age + FC, + data = ma, phenotypes = phen, scalar = "FD", + element.subset = 1:3, + var.terms = c("estimate"), var.model = c("adj.r.squared"), + correct.p.value.terms = c("none"), correct.p.value.model = c("none"), + num.subj.lthr.abs = 0, n_cores = 1, pbar = FALSE, verbose = FALSE + ) + expect_equal(nrow(out), 3L) + + # same fit computed by hand from the in-memory matrices + expected <- vapply(1:3, function(i) { + d <- data.frame(FD = fd[i, ], age = phen$age, FC = fc[i, ]) + unname(stats::coef(stats::lm(FD ~ age + FC, data = d))[["FC"]]) + }, numeric(1)) + expect_equal(out$FC.estimate, expected) +}) + +test_that("chunked TileDB reads match unchunked reads and survive forking", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + old_block <- getOption("ModelArray.tiledb_read_block_size") + on.exit(options(ModelArray.tiledb_read_block_size = old_block), add = TRUE) + + src <- paste0("sub-", 1:6) + set.seed(42) + values <- matrix(stats::rnorm(60), nrow = 10, ncol = 6) + write_tiledb_scalar(store, "FD", values, src) + ma <- ModelArray(store, scalar_types = "FD") + phen <- data.frame(source_file = src, age = seq(10, 60, by = 10)) + + run <- function(element_subset, n_cores) { + ModelArray.lm( + FD ~ age, + data = ma, phenotypes = phen, scalar = "FD", + element.subset = element_subset, + var.terms = c("estimate"), var.model = c("adj.r.squared"), + correct.p.value.terms = c("none"), correct.p.value.model = c("none"), + num.subj.lthr.abs = 0, n_cores = n_cores, pbar = FALSE, verbose = FALSE + ) + } + + options(ModelArray.tiledb_read_block_size = 0L) + unchunked <- run(1:10, 1) + + # block size 3 over 10 elements forces four real chunked reads + options(ModelArray.tiledb_read_block_size = 3L) + expect_equal(.scalar_read_block_size(.build_lm_context(FD ~ age, ma, phen, scalar = "FD")), 3L) + expect_equal(run(1:10, 1), unchunked) + + # out-of-order subset exercises the non-contiguous cache position lookup + shuffled <- as.integer(c(7, 2, 9, 1, 4, 10, 3, 8, 5, 6)) + expect_equal(run(shuffled, 1), unchunked[shuffled, ], ignore_attr = TRUE) + + # forked workers must see the parent's materialized block + skip_on_os("windows") + expect_equal(run(1:10, 2), unchunked) +}) + +test_that("mergeModelArrays records a backend per scalar across storage types", { + skip_if_no_tiledb() + + store <- tempfile(fileext = ".tdb") + h5_file <- tempfile(fileext = ".h5") + on.exit(unlink(store, recursive = TRUE), add = TRUE) + on.exit(unlink(h5_file), add = TRUE) + + src <- paste0("sub-", 1:4) + fd <- matrix(seq_len(12), nrow = 3, ncol = 4) + fc <- matrix(seq_len(12) * 2, nrow = 3, ncol = 4) + write_tiledb_scalar(store, "FD", fd, src) + + h5 <- hdf5r::H5File$new(h5_file, mode = "w") + fc_grp <- h5$create_group("scalars")$create_group("FC") + fc_grp[["values"]] <- fc + fc_grp[["column_names"]] <- src + h5$close_all() + + ma_tiledb <- ModelArray(store, scalar_types = "FD") + ma_hdf5 <- ModelArray(h5_file, scalar_types = "FC") + expect_identical(ma_tiledb@backend, "tiledb") + expect_identical(ma_hdf5@backend, "hdf5") + + phen_tiledb <- data.frame( + subject_id = paste0("subj_", 1:4), + source_file = src, + age = c(10, 20, 30, 40), + stringsAsFactors = FALSE + ) + phen_hdf5 <- data.frame( + subject_id = paste0("subj_", 1:4), + source_file = src, + stringsAsFactors = FALSE + ) + merged <- mergeModelArrays( + list(ma_tiledb, ma_hdf5), + list(phen_tiledb, phen_hdf5), + merge_on = "subject_id" + ) + + expect_identical(merged$data@backend[["FD"]], "tiledb") + expect_identical(merged$data@backend[["FC"]], "hdf5") + # per-scalar lookup resolves through the named path vector + expect_identical(.modelarray_backend(merged$data, merged$data@path[["FD"]]), "tiledb") + expect_identical(.modelarray_backend(merged$data, merged$data@path[["FC"]]), "hdf5") + # a mixed merge still counts as TileDB-attached, so row caching engages + ctx <- .build_lm_context(FD ~ age + FC, merged$data, merged$phenotypes, scalar = "FD") + expect_true(.has_tiledb_attached_scalar(ctx)) +}) + +test_that("auto-detection does not mistake a plain directory for a TileDB store", { + plain_dir <- tempfile() + on.exit(unlink(plain_dir, recursive = TRUE), add = TRUE) + dir.create(plain_dir) + + expect_identical(.resolve_storage_backend(plain_dir, "auto"), "hdf5") + expect_identical(.resolve_storage_backend(file.path(plain_dir, "out.h5"), "auto"), "hdf5") + + # a .tdb suffix declares intent even before the store exists + expect_identical(.resolve_storage_backend(tempfile(fileext = ".tdb"), "auto"), "tiledb") + # ...as does a directory actually shaped like a store + dir.create(file.path(plain_dir, "scalars")) + expect_identical(.resolve_storage_backend(plain_dir, "auto"), "tiledb") + # an explicit choice always wins + expect_identical(.resolve_storage_backend(plain_dir, "hdf5"), "hdf5") +}) diff --git a/tests/testthat/test-wrap_write_scalars.R b/tests/testthat/test-wrap_write_scalars.R index 4b39fbb..8f840bc 100644 --- a/tests/testthat/test-wrap_write_scalars.R +++ b/tests/testthat/test-wrap_write_scalars.R @@ -84,4 +84,22 @@ test_that("ModelArray.wrap write-scalar mode validates column metadata", { ), "length\\(write_scalar_column_names\\) must equal number of selected write_scalar_columns" ) + + expect_error( + ModelArray.wrap( + FUN = simple_fun, + data = modelarray, + phenotypes = phen, + scalar = "FD", + element.subset = as.integer(c(1, 2)), + n_cores = 1, + pbar = FALSE, + verbose = FALSE, + num.subj.lthr.abs = 0, + write_scalar_name = "FD_harmonized", + write_scalar_file = h5_out, + write_scalar_compression_level = 10L + ), + "write_scalar_compression_level must be a single integer between 0 and 9" + ) }) diff --git a/vignettes/container.Rmd b/vignettes/container.Rmd index ade644a..79d7ba6 100644 --- a/vignettes/container.Rmd +++ b/vignettes/container.Rmd @@ -10,13 +10,13 @@ vignette: > ## Preface Target audience: Any user who wants to use a container to run `ModelArray` and/or `ModelArrayIO` commands -What is covered on this page? How to use the container image of `ModelArray + ModelArrayIO`. +What is covered on this page? How to use the container image of `ModelArray + ModelArrayIO`, including the optional TileDB backend. What is not covered on this page? Step by step details on how to run `ModelArray` and `ModelArrayIO` commands. For a full walkthrough on how to use our software, please refer to `vignette("walkthrough")`. We highly suggest reviewing that page if you're new to `ModelArray`. ## Introduction -Besides running `ModelArray` and `ModelArrayIO` commands on a local computer, users also have an option to run them on High Performance Computing (HPC) clusters. When using HPC clusters, users may not have full privilege to install all the dependent packages. Therefore, we provide the option to download our software as a container image that includes `ModelArray + ModelArrayIO`. +Besides running `ModelArray` and `ModelArrayIO` commands on a local computer, users also have an option to run them on High Performance Computing (HPC) clusters. When using HPC clusters, users may not have full privilege to install all the dependent packages. Therefore, we provide the option to download our software as a container image that includes `ModelArray`, `ModelArrayIO`, and the R packages needed for TileDB-backed ModelArray stores. Our container image is publicly available at [Docker Hub](https://hub.docker.com/r/pennlinc/modelarray_confixel). Although HPC clusters do not usually support `docker` commands, this container image can be run with `singularity` commands, which are usually available on HPC clusters (see below for details). @@ -40,6 +40,14 @@ singularity run --cleanenv \ You should now be inside an R environment with the capability to load the `ModelArray` R package as you would in RStudio: `library(ModelArray)`. +You can also check that the optional TileDB backend is available: + +```{.console} +singularity run --cleanenv \ + modelarray_confixel_.sif \ + Rscript -e 'library(ModelArray); stopifnot(requireNamespace("tiledb", quietly = TRUE)); stopifnot(requireNamespace("TileDBArray", quietly = TRUE))' +``` + This is also the way to interactively use this container image to run R, but please make sure you mount necessary directories with `-B` - see below for explanations. ### How to run the container? @@ -55,7 +63,18 @@ Here: `/path/to/singularity/image/modelarray_confixel_.sif` is the full path to the singularity image (`.sif`) you pulled; -`ModelArray` and `ModelArrayIO` are both included in this container. To run conversion tools, simply replace `` with a `ModelArrayIO` command (e.g., `confixel`). +`ModelArray` and `ModelArrayIO` are both included in this container. To convert data, replace `` with the current `ModelArrayIO` CLI, `modelarrayio to-modelarray`. For example, to create a TileDB-backed store: + +```{.console} +singularity run --cleanenv -B /directory/of/your/data \ + /path/to/singularity/image/modelarray_confixel_.sif \ + modelarrayio to-modelarray \ + --backend tiledb \ + --cohort-file /directory/of/your/data/cohort.csv \ + --output /directory/of/your/data/modelarray.tdb +``` + +For HDF5 output, omit `--backend tiledb` or set `--backend hdf5` and use an `.h5` output path. For MIF/fixel data, also provide `--index-file` and `--directions-file`; for NIfTI data, provide `--mask`. To run `ModelArray`, you may first save the R commands into an R script, then replace `` with `Rscript /path/to/your/Rscript`. Make sure this Rscript has been mounted (`-B`), too, e.g., If it's not in the directory `/directory/of/your/data`, please do: `-B /directory/of/your/data,/directory/of/Rscript`. See [here](https://singularity-userdoc.readthedocs.io/en/latest/bind_paths_and_mounts.html) for more on how to mount multiple directories. diff --git a/vignettes/doc_for_developer.Rmd b/vignettes/doc_for_developer.Rmd index df2a9d3..d3d6deb 100644 --- a/vignettes/doc_for_developer.Rmd +++ b/vignettes/doc_for_developer.Rmd @@ -52,9 +52,10 @@ We first define a class called `ModelArray` using `setClass()`. An `ModelArray` * sources: source filenames (e.g., a list of .mif files of fixel-wise data) * scalars: scalar matrix (or matrices) * results: statistical result matrices (if any) -* path: the path to the h5 file on disk +* path: the path to the HDF5 file or TileDB store +* backend: the resolved storage backend (`"hdf5"` or `"tiledb"`) -The key feature of an `ModelArray` object is memory efficient. This is because the entire dataset in HDF5 (.h5) file was not loaded into the memory; only minimal data was loaded. To achieve this, we first need `ModelArraySeed()`, which utilizes `HDF5Array::HDF5ArraySeed()`, acting as a pointer to the .hdf5 file on disk. To make the arrays in the `ModelArray` object look more like "real" arrays, e.g. common array operations such as indexing and transposing can be applied, we utilize package `DelayedArray` to wraps the data in on-disk HDF5 file into a `DelayedArray` object. Finally a `ModelArray` class is defined by integrating above slots together. This is done in the `ModelArray()` function. +The key feature of a `ModelArray` object is memory efficiency: the full dataset is not loaded into memory. `ModelArraySeed()` creates an `HDF5Array::HDF5ArraySeed()` or `TileDBArray::TileDBArraySeed()` that points to the on-disk data. The `DelayedArray` package wraps that seed so common array operations, such as indexing and transposing, remain available. The `ModelArray()` function combines these delayed arrays with the other slots above. Above functions and setups can be found in script `R/ModelArray_Constructor.R` @@ -174,5 +175,3 @@ For more details on how to write unit tests + test out, please check out `Testin The `Dockerfile` for the Docker image can be found in the root folder of the GitHub repository of `ModelArray`. CircleCI will automatically build the Docker image and push to Docker Hub when there is a GitHub commit or a merge to main branch, or when a version is tagged. See the `.circleci/config.yml` file for more details on these setups. - - diff --git a/vignettes/element-splitting.Rmd b/vignettes/element-splitting.Rmd index da35f42..af9b891 100644 --- a/vignettes/element-splitting.Rmd +++ b/vignettes/element-splitting.Rmd @@ -7,6 +7,21 @@ vignette: > %\VignetteEncoding{UTF-8} --- + + ## Why split elements across jobs? ModelArray supports within-process parallelism via the `n_cores` parameter - @@ -39,6 +54,31 @@ for (i in seq_along(starts)) { } ``` +
+TileDB + +Use the TileDB store path when constructing the `ModelArray`. The split range +logic is otherwise unchanged, but the full TileDB version is: + +```{r split_ranges_tiledb, eval=FALSE} +library(ModelArray) + +modelarray <- ModelArray("data.tdb", scalar_types = c("FDC")) +n_total <- numElementsTotal(modelarray, "FDC") + +# Split into chunks of 100,000 elements +chunk_size <- 100000 +starts <- seq(1, n_total, by = chunk_size) +ends <- pmin(starts + chunk_size - 1, n_total) + +# Each job gets one of these ranges +for (i in seq_along(starts)) { + cat(sprintf("Job %d: element.subset = %d:%d\n", i, starts[i], ends[i])) +} +``` + +
+ ```{.console} Job 1: element.subset = 1:100000 Job 2: element.subset = 100001:200000 @@ -84,6 +124,44 @@ chunk_h5 <- sprintf("result_chunk_%d.h5", job_index) writeResults(chunk_h5, df.output = result, analysis_name = "results_lm_chunk") ``` +
+TileDB + +For a TileDB-backed analysis, load `data.tdb` in each job. Save each partial +result to a separate `.rds` file, then combine those files in a single final +write step: + +```{r job_script_tiledb, eval=FALSE} +# run_chunk.R +# Usage: Rscript run_chunk.R + +args <- commandArgs(trailingOnly = TRUE) +job_index <- as.integer(args[1]) + +library(ModelArray) + +modelarray <- ModelArray("data.tdb", scalar_types = c("FDC")) +phenotypes <- read.csv("cohort.csv") + +chunk_size <- 100000 +n_total <- numElementsTotal(modelarray, "FDC") +start <- (job_index - 1) * chunk_size + 1 +end <- min(job_index * chunk_size, n_total) + +result <- ModelArray.lm( + FDC ~ Age + sex + motion, + modelarray, phenotypes, "FDC", + element.subset = start:end, + n_cores = 4 +) + +# Safe because each job writes to a different file +chunk_rds <- sprintf("result_chunk_%d.rds", job_index) +saveRDS(result, chunk_rds) +``` + +
+ On an HPC cluster with a job array (e.g., SLURM): ```{.console} @@ -136,4 +214,47 @@ cat("Total rows:", nrow(full_result), "\n") writeResults("data.h5", df.output = full_result, analysis_name = "results_lm") ``` -**Important**: only one process should write to the HDF5 file at a time. HDF5 does not support concurrent writes from multiple processes — see `vignette("hdf5-large-analyses")` for details on why. +**Important**: only one process should write to the HDF5 file at a time. HDF5 does not support concurrent writes from multiple processes — see `vignette("large-scale-analyses")` for details on why. + +
+TileDB + +After all jobs complete, read each chunk `.rds`, concatenate rows, then write +one final result to `data.tdb`: + +```{r combine_results_tiledb, eval=FALSE} +library(ModelArray) + +chunk_files <- Sys.glob("result_chunk_*.rds") +if (length(chunk_files) == 0) { + stop("No chunk RDS files found (expected files like result_chunk_1.rds).") +} +chunk_ids <- as.integer(sub("result_chunk_([0-9]+)\\.rds$", "\\1", basename(chunk_files))) +if (anyNA(chunk_ids)) { + stop("Could not parse chunk index from one or more chunk filenames.") +} +chunk_files <- chunk_files[order(chunk_ids)] + +results_list <- lapply(chunk_files, readRDS) + +if (!all(vapply(results_list[-1], function(x) identical(colnames(x), colnames(results_list[[1]])), logical(1)))) { + stop("Column names differ across chunks; cannot concatenate safely.") +} + +full_result <- do.call(rbind, results_list) + +# Verify we have all elements +cat("Total rows:", nrow(full_result), "\n") + +writeResults( + "data.tdb", + df.output = full_result, + analysis_name = "results_lm", + backend = "tiledb" +) +``` + +Use the same single-writer pattern for TileDB stores: one final process should +create or replace the analysis result in the shared `.tdb` store. + +
diff --git a/vignettes/elements.Rmd b/vignettes/elements.Rmd index b57fd5a..b7a0f5b 100644 --- a/vignettes/elements.Rmd +++ b/vignettes/elements.Rmd @@ -7,6 +7,21 @@ vignette: > %\VignetteEncoding{UTF-8} --- + + ## What is an element? ModelArray performs mass-univariate statistical analysis - @@ -44,6 +59,27 @@ modelarray <- ModelArray("path/to/data.h5", scalar_types = c("FDC")) scalars(modelarray)[["FDC"]] ``` +
+TileDB + +For TileDB-backed data, the same scalar matrix is stored in a directory-backed +`.tdb` store instead of a single `.h5` file: + +```{r scalar_matrix_tiledb, eval=FALSE} +library(ModelArray) + +# Create a ModelArray object from a TileDB store +modelarray <- ModelArray("path/to/data.tdb", scalar_types = c("FDC")) + +# View the scalar matrix +scalars(modelarray)[["FDC"]] +``` + +ModelArray auto-detects paths ending in `.tdb`; you can also pass +`backend = "tiledb"` explicitly. + +
+ ```{.console} <602229 x 100> matrix of class DelayedMatrix and type "double": FDC/sub-6fee490.mif FDC/sub-647f86c.mif ... FDC/sub-063fd82.mif @@ -54,7 +90,16 @@ scalars(modelarray)[["FDC"]] ``` This is a `DelayedMatrix`, meaning the data lives on disk in the HDF5 file and is only read into memory when accessed. -This is what makes ModelArray memory-efficient — see `vignette("hdf5-large-analyses")` for details. +This is what makes ModelArray memory-efficient — see `vignette("large-scale-analyses")` for details. + +
+TileDB + +TileDB-backed scalar matrices are also exposed through the `DelayedArray` +interface. Values are read from the `.tdb` store only when requested, so the +same accessor and modelling code works for both backends. + +
## Element IDs diff --git a/vignettes/exploring-h5.Rmd b/vignettes/exploring-modelarray-data.Rmd similarity index 61% rename from vignettes/exploring-h5.Rmd rename to vignettes/exploring-modelarray-data.Rmd index 295a4b9..914504a 100644 --- a/vignettes/exploring-h5.Rmd +++ b/vignettes/exploring-modelarray-data.Rmd @@ -1,26 +1,72 @@ --- -title: "Exploring HDF5 Data with Convenience Functions" +title: "Exploring ModelArray Data with Convenience Functions" output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{Exploring HDF5 Data with Convenience Functions} + %\VignetteIndexEntry{Exploring ModelArray Data with Convenience Functions} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ModelArray provides several accessor functions for inspecting your data without writing analysis code. -This vignette walks through each one. +This vignette walks through each one for either HDF5 files or TileDB stores. +The main examples use HDF5, with TileDB-specific setup shown in blue toggles where the code differs. + + ```{r setup, eval=FALSE} library(ModelArray) -h5_path <- "~/Desktop/myProject/demo_FDC_n100.h5" -modelarray <- ModelArray(h5_path, +modelarray_path <- "~/Desktop/myProject/demo_FDC_n100.h5" +csv_path <- "~/Desktop/myProject/cohort_FDC_n100.csv" + +modelarray <- ModelArray(modelarray_path, scalar_types = c("FDC"), analysis_names = c("results_lm") ) -phenotypes <- read.csv("~/Desktop/myProject/cohort_FDC_n100.csv") +phenotypes <- read.csv(csv_path) ``` +
+TileDB + +For a TileDB-backed store created by `modelarrayio to-modelarray --backend tiledb`, +use the `.tdb` directory path. ModelArray auto-detects `.tdb` paths; you can also +pass `backend = "tiledb"` explicitly. + +```{r setup_tiledb, eval=FALSE} +library(ModelArray) + +modelarray_path <- "~/Desktop/myProject/demo_FDC_n100.tdb" +csv_path <- "~/Desktop/myProject/cohort_FDC_n100.csv" + +modelarray <- ModelArray(modelarray_path, + scalar_types = c("FDC"), + analysis_names = c("results_lm") +) +phenotypes <- read.csv(csv_path) +``` + +TileDB support in R requires the optional `tiledb` and `TileDBArray` packages. +See the optional TileDB installation notes in `vignette("installations")`. + +
+ +The accessor examples below start from `modelarray`, so they are identical for +HDF5 and TileDB once the object has been constructed. + ## `show()` — Object summary Simply typing the object name (or calling `show()`) prints a summary: @@ -33,12 +79,27 @@ modelarray ModelArray located at ~/Desktop/myProject/demo_FDC_n100.h5 Source files: 100 - Scalars: FDC + Scalars: FDC (602229 elements) Analyses: results_lm ``` This tells you at a glance: -the file path, how many sources, which scalars are loaded, and which analyses have been saved. +the storage path, how many sources, which scalars are loaded, and which analyses have been saved. + +
+TileDB + +For a TileDB-backed object, the printed path points to the `.tdb` store: + +```{.console} +ModelArray located at ~/Desktop/myProject/demo_FDC_n100.tdb + + Source files: 100 + Scalars: FDC (602229 elements) + Analyses: results_lm +``` + +
## `sources()` — Source files @@ -64,7 +125,7 @@ Access the full scalar matrix or a specific scalar: # All scalars (returns a named list) scalars(modelarray) -# A specific scalar (returns a DelayedMatrix) +# A specific scalar (returns a DelayedArray-backed matrix) fdc_matrix <- scalars(modelarray)[["FDC"]] dim(fdc_matrix) ``` @@ -89,7 +150,7 @@ subject_values <- as.numeric(fdc_matrix[, 1]) length(subject_values) ``` -Note: because the data is stored on disk as a `DelayedMatrix`, accessing a row or column triggers an HDF5 read. +Note: because the data is stored on disk and exposed through `DelayedArray`, accessing a row or column triggers a backend read from HDF5 or TileDB. Avoid looping over many rows manually: that's what `ModelArray.lm()` and friends are optimized for. ## `results()` — Saved analysis results @@ -101,7 +162,7 @@ If you loaded analysis names when creating the object, you can access them: results(modelarray) # A specific analysis -lm_results <- results(modelarray)[["results_lm"]] +lm_results <- results(modelarray)[["results_lm"]]$results_matrix dim(lm_results) colnames(lm_results) ``` @@ -111,31 +172,51 @@ make sure you specified it in the `analysis_names` argument when calling `ModelA ```{r results_load, eval=FALSE} # This loads results; without analysis_names, the results slot is empty -modelarray <- ModelArray(h5_path, +modelarray <- ModelArray(modelarray_path, scalar_types = "FDC", analysis_names = c("results_lm") ) ``` -## Quick inspection with `h5summary()` +## Quick inspection with `ModelArraySummary()` -Before constructing a full ModelArray, you can inspect an h5 file's structure: +Before constructing a full ModelArray, you can inspect the storage structure: ```{r h5summary, eval=FALSE} -h5summary(h5_path) +ModelArraySummary(modelarray_path) +``` + +```{.console} +H5 file: ~/Desktop/myProject/demo_FDC_n100.h5 + +Scalars: + FDC: 602229 elements x 100 input files + +Analyses: results_lm +``` + +
+TileDB + +Use the TileDB store path the same way. The explicit backend argument is optional +for `.tdb` paths: + +```{r tiledbsummary, eval=FALSE} +ModelArraySummary("~/Desktop/myProject/demo_FDC_n100.tdb", backend = "tiledb") ``` ```{.console} -ModelArray HDF5 summary: ~/Desktop/myProject/demo_FDC_n100.h5 +TileDB store: ~/Desktop/myProject/demo_FDC_n100.tdb - Scalars: - name nElements nInputFiles - 1: FDC 602229 100 +Scalars: + FDC: 602229 elements x 100 input files - Analyses: results_lm +Analyses: results_lm ``` -This is lightweight — it reads the h5 metadata without loading any data into memory. +
+ +This is lightweight — it reads storage metadata without loading any data into memory. ## Dimension accessors @@ -159,7 +240,7 @@ The older `numElementsTotal()` function still works but `nElements()` and `nInpu ## `elementMetadata()` — Spatial metadata -Some h5 files contain per-element metadata (e.g., greyordinate labels for cifti data, +Some ModelArray stores contain per-element metadata (e.g., greyordinate labels for cifti data, fixel directions, or voxel coordinates). Access it with: ```{r element_metadata, eval=FALSE} @@ -168,7 +249,7 @@ dim(em) # e.g. 602229 x 3 head(em) ``` -Returns `NULL` if the h5 file does not contain element metadata. +Returns `NULL` if the store does not contain element metadata. ## `numElementsTotal()` — Element count (legacy) @@ -230,7 +311,7 @@ sum(!is.na(element_values)) After running a model, you can query the saved results to find the most significant elements: ```{r strongest_effect, eval=FALSE} -lm_results <- results(modelarray)[["results_lm"]] +lm_results <- results(modelarray)[["results_lm"]]$results_matrix # Find the column index for Age p-value (FDR-corrected) col_idx <- which(colnames(lm_results) == "Age.p.value.fdr") diff --git a/vignettes/hdf5-format.Rmd b/vignettes/hdf5-format.Rmd deleted file mode 100644 index 70efea4..0000000 --- a/vignettes/hdf5-format.Rmd +++ /dev/null @@ -1,134 +0,0 @@ ---- -title: "Understanding the HDF5 File" -output: rmarkdown::html_vignette -vignette: > - %\VignetteIndexEntry{Understanding the HDF5 File} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -## What is HDF5? - -HDF5 (Hierarchical Data Format version 5) is a file format designed for storing large, structured datasets. -Think of it as a filesystem within a file: - -- **Groups** act like directories, organizing data hierarchically -- **Datasets** act like files, holding arrays of data -- **Attributes** act like metadata, attached to groups or datasets - -HDF5 files use the `.h5` extension and are widely used in scientific computing because they support on-disk access -(reading data without loading the entire file into memory), compression, and efficient chunked storage. - -## The ModelArray HDF5 layout - -ModelArray expects a specific structure inside the HDF5 file. -This structure is created by companion tools — -[ModelArrayIO](https://github.com/PennLINC/ModelArrayIO) -for fixel data, voxel data, and CIFTI data. - -The layout follows this pattern: - -``` -/ -├── scalars/ -│ └── / # e.g., "FDC", "FA", "thickness" -│ └── values # Dataset: elements × source files matrix -│ └── column_names # Dataset (text): source filenames -└── results/ # Created by writeResults() - └── / # e.g., "results_lm", "results_gam" - ├── results_matrix # Dataset: elements × statistics matrix - └── column_names # Dataset (text): statistic column names -``` - -### The scalars group - -`/scalars//values` is the main data matrix: - -- **Rows** = elements (fixels, voxels, or greyordinates) -- **Columns** = source files - -`/scalars//column_names` stores the source filenames that correspond to each column, -so you can trace every column back to its subject. - -In modern `ModelArrayIO`, column names are written to a dedicated text dataset because long name vectors can exceed practical HDF5 attribute limits. - -### The results group - -`/results//results_matrix` is created when you call `writeResults()` to save statistical outputs. -Each analysis gets its own subgroup, so you can store multiple analyses -(e.g., `results_lm` and `results_gam`) in the same file. - -The corresponding statistic names are stored in `/results//column_names` as a text dataset. - -### Modern vs legacy column-name storage - -- **Preferred (current)**: text dataset `column_names` -- **Legacy (older files)**: attributes on `values`/`results_matrix` (for example `column_names` or `colnames`) - -If both are present, prefer the `column_names` dataset as the canonical source and treat attributes as backward-compatibility metadata. - -## How ModelArray reads HDF5 - -When you create a ModelArray object, the scalar data is **not** loaded into memory. -Instead, ModelArray uses the [DelayedArray](https://bioconductor.org/packages/DelayedArray/) framework to create a lazy reference to the on-disk data: - -```{r delayed_array, eval=FALSE} -library(ModelArray) - -modelarray <- ModelArray("data.h5", scalar_types = c("FDC")) -scalars(modelarray)[["FDC"]] -``` - -```{.console} -<602229 x 100> matrix of class DelayedMatrix and type "double": -``` - -The `DelayedMatrix` holds a pointer to the HDF5 file. Data is only read from disk when you actually access specific rows or columns. -This is why ModelArray can handle datasets with hundreds of thousands of elements and thousands of subjects without running out of memory. - -During model fitting, ModelArray reads **one element (row) at a time** — pulling a single row of ~N subject values, fitting the model, and moving on. At no point is the full matrix loaded into RAM. - -## Inspecting an HDF5 file - -You can explore the structure of any HDF5 file using `rhdf5::h5ls()`: - -```{r inspect_h5, eval=FALSE} -rhdf5::h5ls("data.h5") -``` - -```{.console} - group name otype dclass dim -0 / analysis_configs H5I_GROUP -1 /analysis_configs results_lm H5I_GROUP -2 / results H5I_GROUP -3 /results results_lm H5I_GROUP -4 /results/results_lm results_matrix H5I_DATASET FLOAT 602229 x 17 -5 / scalars H5I_GROUP -6 /scalars FDC H5I_GROUP -7 /scalars/FDC values H5I_DATASET FLOAT 602229 x 100 -``` - -You can also read specific pieces of data directly: - -```{r read_h5, eval=FALSE} -# Preferred: read names from the text dataset -rhdf5::h5read("data.h5", "scalars/FDC/column_names") - -# Legacy fallback for older files -rhdf5::h5readAttributes("data.h5", "scalars/FDC/values")$column_names - -# Read a small slice of the data matrix -rhdf5::h5read("data.h5", "scalars/FDC/values", index = list(1:5, 1:3)) -``` - -## Creating HDF5 files - -ModelArray does not create HDF5 files from raw imaging data — that's the job of the companion conversion tools: - -| Data type | Command | -|:----------|:--------| -| Fixel (`.mif`) | `confixel` | -| Voxel (`.nii.gz`) | `convoxel` | -| Surface (`.dscalar.nii`) | `concifti` | - -Each command reads the source imaging files listed in a cohort CSV and writes them into the HDF5 layout described above. See the [ModelArrayIO documentation](https://github.com/PennLINC/ModelArrayIO) for detailed usage instructions. diff --git a/vignettes/hdf5-large-analyses.Rmd b/vignettes/hdf5-large-analyses.Rmd deleted file mode 100644 index 7d211ac..0000000 --- a/vignettes/hdf5-large-analyses.Rmd +++ /dev/null @@ -1,99 +0,0 @@ ---- -title: "Using HDF5 for Large-Scale Analyses" -output: rmarkdown::html_vignette -vignette: > - %\VignetteIndexEntry{Using HDF5 for Large-Scale Analyses} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -## Why HDF5? - -A typical fixel-based analysis might involve 600,000 fixels across 1,000 subjects. Stored naively as a double-precision matrix, that's about 4.5 GB — too large to hold in memory on most systems, let alone alongside the R session overhead and model fitting. - -HDF5 solves this by keeping data on disk and providing efficient random access to slices of the data. ModelArray leverages this to fit models element-by-element without ever loading the full matrix into RAM. - -## Chunking - -HDF5 files store data in **chunks** — rectangular blocks of the array that are read and written as a unit. When you access a single row (element), HDF5 only reads the chunk(s) that contain that row, not the entire dataset. - -The chunk layout matters for performance. Because ModelArray reads data **row by row** (one element at a time), chunks that span a modest number of rows and all columns work best. This is the default behavior of the `ModelArrayIO` conversion commands (`confixel`, `convoxel`, `concifti`). - -When creating HDF5 files with the `confixel` command from `ModelArrayIO`, you can control the chunk size with the `--chunkmb` flag: - -```{.console} -$ confixel \ - --index-file FDC/index.mif \ - --directions-file FDC/directions.mif \ - --cohort-file cohort.csv \ - --relative-root /path/to/project \ - --output-hdf5 data.h5 \ - --chunkmb 4 # chunk size in MB (default) -``` - -Smaller chunks use less memory per read but may require more disk seeks. The default of 4 MB is a reasonable trade-off for most datasets. - -## Compression - -HDF5 supports transparent compression of chunks using gzip. When enabled, data is compressed on disk and decompressed on-the-fly when read. This typically reduces file sizes by 30–60% with minimal impact on read speed. - -`ModelArrayIO` enables gzip compression by default. You can see this reflected in the output filenames: - -```{.console} -study-HBN_compression-gzip_chunkmb-4_thickness.h5 -``` - -The compression is transparent to ModelArray — it reads the data the same way regardless of whether it's compressed. - -## Memory efficiency during model fitting - -When you call `ModelArray.lm()`, `ModelArray.gam()`, or `ModelArray.wrap()`, here is what happens for each element: - -1. **Read one row** from the scalar matrix (the values for all subjects at this element) -2. **Combine** with the phenotypes data frame to create a per-element data frame -3. **Fit the model** (lm, gam, or your custom function) -4. **Extract statistics** into a single row of the output data frame -5. **Discard** the per-element data and move to the next element - -At no point is more than one row of the scalar matrix in memory. This is why ModelArray can handle arbitrarily large datasets — the memory footprint is determined by the number of subjects (columns), not the number of elements (rows). - -## Parallelism within a single process - -ModelArray supports parallel processing via the `n_cores` parameter: - -```{r parallel, eval=FALSE} -result <- ModelArray.lm(FDC ~ Age + sex, modelarray, phenotypes, "FDC", - n_cores = 4 -) -``` - -This uses `parallel::mclapply()` (fork-based parallelism on Linux/macOS) to process multiple elements simultaneously. Each forked worker inherits the same file handle and reads different rows from the HDF5 file. Because HDF5 supports concurrent **reads** from the same file, this works safely and scales well. - -## Why concurrent writes don't work - -While multiple processes can safely **read** from the same HDF5 file, they **cannot safely write** to it simultaneously. The standard HDF5 library (without special MPI-IO compilation) does not support concurrent writes — doing so can corrupt the file. - -This means you should **never** have multiple R processes calling `writeResults()` on the same HDF5 file at the same time. If you split your analysis across HPC jobs (see `vignette("element-splitting")`), each job should save its partial results to a separate `.rds` file, then a single final step combines them and writes to the HDF5 file. - -```{r safe_write, eval=FALSE} -# WRONG: multiple jobs writing to the same H5 simultaneously -# writeResults("data.h5", df.output = my_partial_result, ...) # Don't do this in parallel! - -# RIGHT: save partial results, then combine in one process -saveRDS(my_partial_result, sprintf("result_chunk_%d.rds", job_id)) - -# ... after all jobs finish, in a single process: -chunks <- lapply(Sys.glob("result_chunk_*.rds"), readRDS) -full_result <- do.call(rbind, chunks) -writeResults("data.h5", df.output = full_result, analysis_name = "results_lm") -``` - -## Backing up your HDF5 file - -Because `writeResults()` modifies the HDF5 file in place, it's good practice to back up your file before writing results: - -```{.console} -$ cp data.h5 data_backup.h5 -``` - -If something goes wrong during a write (e.g., the process is killed mid-write), the HDF5 file may be left in an inconsistent state. Having a backup lets you recover without re-running the conversion step. diff --git a/vignettes/installations.Rmd b/vignettes/installations.Rmd index 9f99d9b..71bbf28 100644 --- a/vignettes/installations.Rmd +++ b/vignettes/installations.Rmd @@ -20,19 +20,19 @@ We will first set up the conda environment, then install some dependent packages We first create a conda environment `modelarray` for installing the companion software `ModelArrayIO` etc. We'll install python version 3.9: ``` {.console} -foo@bar:~$ conda create --name modelarray python=3.9 +foo@bar:~$ conda create --name modelarray python=3.11 foo@bar:~$ conda activate modelarray ``` ### Install MRtrix (Only required for fixel-wise data) -When converting fixel-wise data (`.mif`), the `confixel` command from `ModelArrayIO` uses `mrconvert` from MRtrix, so please make sure MRtrix has been installed. It can either be installed via `conda` in this conda environment we just created, or be compiled from source. See [MRtrix's webpage](https://www.mrtrix.org/download/) for more. Type `mrview` in the terminal to check whether MRtrix installation is successful. +When working with fixel-wise data (`.mif`), install MRtrix so you can inspect and visualize MIF inputs and outputs. It can either be installed via `conda` in this conda environment we just created, or be compiled from source. See [MRtrix's webpage](https://www.mrtrix.org/download/) for more. Type `mrview` in the terminal to check whether MRtrix installation is successful. If your input data is voxel-wise data, you can skip this step. ### Install HDF5 libraries in the system -Because ModelArray works with the Hierarchical Data Format 5 (HDF5) file format, we need to make sure necessary libraries of HDF5 are installed in the system. +ModelArray always supports data backed by the Hierarchical Data Format 5 (HDF5) file format, so we need to make sure the necessary HDF5 libraries are installed in the system. TileDB support is optional and requires additional R packages; see the optional TileDB section below if you plan to use TileDB-backed stores. #### On a Linux Ubuntu system @@ -60,7 +60,7 @@ foo@bar:~$ brew install hdf5 For details you may refer to the webpage [here](https://formulae.brew.sh/formula/hdf5) ### Install ModelArrayIO python package from GitHub -[ModelArrayIO](https://github.com/PennLINC/ModelArrayIO) provides file format conversion commands (`confixel`, `convoxel`, `concifti`) for fixel-wise data (`.mif`), voxel-wise data (NIfTI), and CIFTI data. Follow the commands below to install it from GitHub: +[ModelArrayIO](https://github.com/PennLINC/ModelArrayIO) provides the `modelarrayio` command-line tool for converting fixel-wise data (`.mif`), voxel-wise data (NIfTI), and CIFTI data into ModelArray storage, and for exporting saved ModelArray results back to neuroimaging formats. Follow the commands below to install it from GitHub: ``` {.console} # We first activate the conda environment we just created: @@ -77,6 +77,12 @@ foo@bar:myProject$ cd .. foo@bar:myProject$ rm -r ModelArrayIO ``` +After installation, conversion is done with `modelarrayio to-modelarray`. The command auto-detects the imaging modality from the source files in the cohort CSV. For example: + +``` {.console} +foo@bar:myProject$ modelarrayio to-modelarray --help +``` + ### Install R @@ -99,6 +105,33 @@ library(devtools) devtools::install_github("PennLINC/ModelArray") ``` +### (Optional) Install TileDB backend dependencies + +ModelArray can also read and write TileDB-backed stores. This is optional: you only need these packages if you plan to use `backend = "tiledb"` or paths ending in `.tdb`. + +```{r install-tiledb, eval=FALSE} +install.packages(c("jsonlite", "tiledb")) + +if (!requireNamespace("BiocManager", quietly = TRUE)) { + install.packages("BiocManager") +} +BiocManager::install("TileDBArray") +``` + +You can verify TileDB support with: + +```{r check-tiledb, eval=FALSE} +library(ModelArray) + +ModelArraySummary("path/to/store.tdb", backend = "tiledb") +``` + +If installing the `tiledb` R package fails, install the TileDB system library first. On macOS: + +``` {.console} +foo@bar:~$ brew install tiledb +``` + Now, ModelArray is ready to use: ```{r} library(ModelArray) diff --git a/vignettes/large-scale-analyses.Rmd b/vignettes/large-scale-analyses.Rmd new file mode 100644 index 0000000..c3470cb --- /dev/null +++ b/vignettes/large-scale-analyses.Rmd @@ -0,0 +1,244 @@ +--- +title: "Large-Scale Analyses" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Large-Scale Analyses} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + + + +## Why HDF5? + +A typical fixel-based analysis might involve 600,000 fixels across 1,000 subjects. Stored naively as a double-precision matrix, that's about 4.5 GB — too large to hold in memory on most systems, let alone alongside the R session overhead and model fitting. + +HDF5 solves this by keeping data on disk and providing efficient random access to slices of the data. ModelArray leverages this to fit models element-by-element without ever loading the full matrix into RAM. + +
+TileDB + +TileDB solves the same on-disk storage problem with a directory-backed `.tdb` +store instead of a single `.h5` file. ModelArray uses the same delayed access +pattern for TileDB-backed stores, so the modelling workflow stays the same once +the store is loaded. + +
+ +## Chunking + +HDF5 files store data in **chunks** — rectangular blocks of the array that are read and written as a unit. When you access a single row (element), HDF5 only reads the chunk(s) that contain that row, not the entire dataset. + +The chunk layout matters for performance. Because ModelArray reads data **row by row** (one element at a time), chunks that span a modest number of rows and all columns work best. This is the default behavior of the `modelarrayio to-modelarray` conversion command. + +When creating HDF5 files with `modelarrayio to-modelarray`, you can control the target chunk size with the `--target-chunk-mb` flag: + +```{.console} +$ modelarrayio to-modelarray \ + --backend hdf5 \ + --index-file FDC/index.mif \ + --directions-file FDC/directions.mif \ + --cohort-file cohort.csv \ + --output data.h5 \ + --target-chunk-mb 4 # target chunk size in MiB +``` + +
+TileDB + +Use `--backend tiledb` and a `.tdb` output path. In TileDB, the +storage-level analogue of an HDF5 chunk is a **tile**. ModelArrayIO uses the +same `--target-chunk-mb` option to choose an approximate TileDB tile size: + +```{.console} +$ modelarrayio to-modelarray \ + --backend tiledb \ + --index-file FDC/index.mif \ + --directions-file FDC/directions.mif \ + --cohort-file cohort.csv \ + --output data.tdb \ + --target-chunk-mb 4 # target TileDB tile size in MiB +``` + +This controls the on-disk TileDB layout. It is separate from ModelArray's +analysis-time TileDB read-block settings, described below under memory +efficiency. + +
+ +For HDF5, smaller chunks use less memory per read but may require more disk seeks. +The default target is 2 MiB; the examples above explicitly override it to 4 MiB. + +## Compression + +HDF5 supports transparent compression of chunks using gzip. When enabled, data is compressed on disk and decompressed on-the-fly when read. This typically reduces file sizes by 30–60% with minimal impact on read speed. + +`ModelArrayIO` enables gzip compression by default. You can see this reflected in the output filenames: + +```{.console} +study-HBN_compression-gzip_chunkmb-4_thickness.h5 +``` + +The compression is transparent to ModelArray — it reads the data the same way regardless of whether it's compressed. + +
+TileDB + +TileDB stores are directories, so the output path ends in `.tdb`: + +```{.console} +study-HBN_compression-gzip_chunkmb-4_thickness.tdb +``` + +ModelArrayIO's compression options map to TileDB compression filters and tile +sizes. As with HDF5, compression is transparent when you load the store in R. + +
+ +## Memory efficiency during model fitting + +When you call `ModelArray.lm()`, `ModelArray.gam()`, or `ModelArray.wrap()`, here is what happens for each element: + +1. **Read one row** from the scalar matrix (the values for all subjects at this element) +2. **Combine** with the phenotypes data frame to create a per-element data frame +3. **Fit the model** (lm, gam, or your custom function) +4. **Extract statistics** into a single row of the output data frame +5. **Discard** the per-element data and move to the next element + +At no point is more than one row of the scalar matrix in memory. This is why ModelArray can handle arbitrarily large datasets — the memory footprint is determined by the number of subjects (columns), not the number of elements (rows). + +
+TileDB + +For TileDB-backed data, ModelArray may hold a bounded block of rows in memory +to avoid slow one-row TileDB reads. The memory footprint is still controlled: it +depends on the number of subjects, the number of attached scalars, and the +configured block size, not the full number of elements in the store. + +The default TileDB read-block target is 512 MiB. Tune it when needed: + +```{r tiledb_block_options, eval=FALSE} +# Target roughly 256 MiB per TileDB read block +options(ModelArray.tiledb_read_block_mb = 256) + +# Or set an explicit number of elements per block +options(ModelArray.tiledb_read_block_size = 5000) + +# Disable TileDB row-block materialization and read row by row +options(ModelArray.tiledb_read_block_size = 0) +``` + +Approximate memory for one TileDB read block is: + +```{.text} +block_size x number_of_subjects x number_of_attached_scalars x 8 bytes +``` + +For example, a block of 5,000 elements across 1,000 subjects for one scalar is +about 40 MB before R object overhead. + +
+ +## Parallelism within a single process + +ModelArray supports parallel processing via the `n_cores` parameter: + +```{r parallel, eval=FALSE} +result <- ModelArray.lm(FDC ~ Age + sex, modelarray, phenotypes, "FDC", + n_cores = 4 +) +``` + +This uses `parallel::mclapply()` (fork-based parallelism on Linux/macOS) to process multiple elements simultaneously. Each forked worker inherits the same file handle and reads different rows from the HDF5 file. Because HDF5 supports concurrent **reads** from the same file, this works safely and scales well. + +
+TileDB + +For TileDB-backed scalars, ModelArray reads each row block in the parent process +before dispatching elements to workers. On Linux/macOS this lets forked workers +reuse the in-memory block through copy-on-write, reducing repeated TileDB reads. + +
+ +## Why concurrent writes don't work + +While multiple processes can safely **read** from the same HDF5 file, they **cannot safely write** to it simultaneously. The standard HDF5 library (without special MPI-IO compilation) does not support concurrent writes — doing so can corrupt the file. + +This means you should **never** have multiple R processes calling `writeResults()` on the same HDF5 file at the same time. If you split your analysis across HPC jobs (see `vignette("element-splitting")`), each job should save its partial results to a separate `.rds` file, then a single final step combines them and writes to the HDF5 file. + +```{r safe_write, eval=FALSE} +# WRONG: multiple jobs writing to the same H5 simultaneously +# writeResults("data.h5", df.output = my_partial_result, ...) # Don't do this in parallel! + +# RIGHT: save partial results, then combine in one process +saveRDS(my_partial_result, sprintf("result_chunk_%d.rds", job_id)) + +# ... after all jobs finish, in a single process: +chunks <- lapply(Sys.glob("result_chunk_*.rds"), readRDS) +full_result <- do.call(rbind, chunks) +writeResults("data.h5", df.output = full_result, analysis_name = "results_lm") +``` + +
+TileDB + +Use the same single-writer pattern for ModelArray writes to TileDB stores. Even +though TileDB has its own concurrency features, `writeResults()` is intended to +create or replace one analysis result from one R process at a time. + +```{r safe_tiledb_write, eval=FALSE} +# WRONG: multiple jobs writing to the same TileDB store simultaneously +# writeResults("data.tdb", df.output = my_partial_result, backend = "tiledb") + +# RIGHT: save partial results, then combine in one process +saveRDS(my_partial_result, sprintf("result_chunk_%d.rds", job_id)) + +# ... after all jobs finish, in a single process: +chunks <- lapply(Sys.glob("result_chunk_*.rds"), readRDS) +full_result <- do.call(rbind, chunks) +writeResults("data.tdb", + df.output = full_result, + analysis_name = "results_lm", + backend = "tiledb" +) +``` + +
+ +## Backing up your HDF5 file + +Because `writeResults()` modifies the HDF5 file in place, it's good practice to back up your file before writing results: + +```{.console} +$ cp data.h5 data_backup.h5 +``` + +If something goes wrong during a write (e.g., the process is killed mid-write), the HDF5 file may be left in an inconsistent state. Having a backup lets you recover without re-running the conversion step. + +
+TileDB + +TileDB stores are directories, so back up the full directory before writing +results: + +```{.console} +$ cp -R data.tdb data_backup.tdb +``` + +If a TileDB write is interrupted, restore from the backup directory before +rerunning the write step. + +
diff --git a/vignettes/modelarray-storage.Rmd b/vignettes/modelarray-storage.Rmd new file mode 100644 index 0000000..e6d9c48 --- /dev/null +++ b/vignettes/modelarray-storage.Rmd @@ -0,0 +1,261 @@ +--- +title: "Understanding ModelArray Storage" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Understanding ModelArray Storage} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + + + +## What is HDF5? + +HDF5 (Hierarchical Data Format version 5) is a file format designed for storing large, structured datasets. +Think of it as a filesystem within a file: + +- **Groups** act like directories, organizing data hierarchically +- **Datasets** act like files, holding arrays of data +- **Attributes** act like metadata, attached to groups or datasets + +HDF5 files use the `.h5` extension and are widely used in scientific computing because they support on-disk access +(reading data without loading the entire file into memory), compression, and efficient chunked storage. + +
+TileDB + +TileDB is an optional backend for storing the same ModelArray content in a +directory-backed `.tdb` store. A TileDB store is organized as groups and arrays +rather than groups and datasets, but it provides the same core features that +ModelArray needs: on-disk access, compression, and tiled array storage. + +
+ +## The ModelArray storage layout + +ModelArray expects a specific structure inside the HDF5 file. +This structure is created by companion tools — +[ModelArrayIO](https://github.com/PennLINC/ModelArrayIO) +for fixel data, voxel data, and CIFTI data. + +The layout follows this pattern: + +``` +/ +├── scalars/ +│ └── / # e.g., "FDC", "FA", "thickness" +│ └── values # Dataset: elements × source files matrix +│ └── column_names # Dataset (text): source filenames +└── results/ # Created by writeResults() + └── / # e.g., "results_lm", "results_gam" + ├── results_matrix # Dataset: elements × statistics matrix + └── column_names # Dataset (text): statistic column names +``` + +
+TileDB + +TileDB-backed stores use the same conceptual layout, represented by TileDB +groups and arrays under a `.tdb` directory: + +```{.text} +data.tdb/ +├── scalars/ +│ └── / # e.g., "FDC", "FA", "thickness" +│ ├── values # TileDB array: elements x source files matrix +│ └── column_names # Optional companion array or metadata +└── results/ # Created by writeResults() + └── / # e.g., "results_lm", "results_gam" + ├── results_matrix # TileDB array: elements x statistics matrix + └── column_names # Optional companion array or metadata +``` + +
+ +### The scalars group + +`/scalars//values` is the main data matrix: + +- **Rows** = elements (fixels, voxels, or greyordinates) +- **Columns** = source files + +`/scalars//column_names` stores the source filenames that correspond to each column, +so you can trace every column back to its subject. + +In modern `ModelArrayIO`, column names are written to a dedicated text dataset because long name vectors can exceed practical HDF5 attribute limits. + +### The results group + +`/results//results_matrix` is created when you call `writeResults()` to save statistical outputs. +Each analysis gets its own subgroup, so you can store multiple analyses +(e.g., `results_lm` and `results_gam`) in the same file. + +The corresponding statistic names are stored in `/results//column_names` as a text dataset. + +### Modern vs legacy column-name storage + +- **Preferred (current)**: text dataset `column_names` +- **Legacy (older files)**: attributes on `values`/`results_matrix` (for example `column_names` or `colnames`) + +If both are present, prefer the `column_names` dataset as the canonical source and treat attributes as backward-compatibility metadata. + +
+TileDB + +TileDB-backed stores use metadata and/or companion arrays for column names. +ModelArray handles the supported TileDB locations internally; most users should +use `sources()`, `scalarNames()`, `analysisNames()`, and `results()` rather than +reading TileDB metadata directly. + +```{r tiledb_column_name_accessors, eval=FALSE} +library(ModelArray) + +modelarray <- ModelArray( + "data.tdb", + scalar_types = c("FDC"), + analysis_names = c("results_lm") +) + +# Scalar names and source filenames +scalarNames(modelarray) +head(sources(modelarray)[["FDC"]]) + +# Saved analysis names and result column names +analysisNames(modelarray) +lm_results <- results(modelarray)[["results_lm"]]$results_matrix +colnames(lm_results) +``` + +
+ +## How ModelArray reads HDF5 + +When you create a ModelArray object, the scalar data is **not** loaded into memory. +Instead, ModelArray uses the [DelayedArray](https://bioconductor.org/packages/DelayedArray/) framework to create a lazy reference to the on-disk data: + +```{r delayed_array, eval=FALSE} +library(ModelArray) + +modelarray <- ModelArray("data.h5", scalar_types = c("FDC")) +scalars(modelarray)[["FDC"]] +``` + +```{.console} +<602229 x 100> matrix of class DelayedMatrix and type "double": +``` + +The `DelayedMatrix` holds a pointer to the HDF5 file. Data is only read from disk when you actually access specific rows or columns. +This is why ModelArray can handle datasets with hundreds of thousands of elements and thousands of subjects without running out of memory. + +During model fitting, ModelArray reads **one element (row) at a time** — pulling a single row of ~N subject values, fitting the model, and moving on. At no point is the full matrix loaded into RAM. + +
+TileDB + +TileDB-backed stores are exposed through the same `DelayedArray` interface via +`TileDBArray`: + +```{r delayed_array_tiledb, eval=FALSE} +library(ModelArray) + +modelarray <- ModelArray("data.tdb", scalar_types = c("FDC")) +scalars(modelarray)[["FDC"]] +``` + +Values are still read lazily from disk. During model fitting, ModelArray may +materialize bounded row blocks for TileDB-backed scalars to avoid slow one-row +TileDB reads; see `vignette("large-scale-analyses")` for details. + +
+ +## Inspecting an HDF5 file + +You can explore the structure of any HDF5 file using `rhdf5::h5ls()`: + +```{r inspect_h5, eval=FALSE} +rhdf5::h5ls("data.h5") +``` + +```{.console} + group name otype dclass dim +0 / analysis_configs H5I_GROUP +1 /analysis_configs results_lm H5I_GROUP +2 / results H5I_GROUP +3 /results results_lm H5I_GROUP +4 /results/results_lm results_matrix H5I_DATASET FLOAT 602229 x 17 +5 / scalars H5I_GROUP +6 /scalars FDC H5I_GROUP +7 /scalars/FDC values H5I_DATASET FLOAT 602229 x 100 +``` + +You can also read specific pieces of data directly: + +```{r read_h5, eval=FALSE} +# Preferred: read names from the text dataset +rhdf5::h5read("data.h5", "scalars/FDC/column_names") + +# Legacy fallback for older files +rhdf5::h5readAttributes("data.h5", "scalars/FDC/values")$column_names + +# Read a small slice of the data matrix +rhdf5::h5read("data.h5", "scalars/FDC/values", index = list(1:5, 1:3)) +``` + +
+TileDB + +Use ModelArray accessors for backend-independent inspection: + +```{r inspect_tiledb, eval=FALSE} +library(ModelArray) + +ModelArraySummary("data.tdb", backend = "tiledb") + +modelarray <- ModelArray("data.tdb", scalar_types = c("FDC")) +sources(modelarray)[["FDC"]] +as.matrix(scalars(modelarray)[["FDC"]][1:5, 1:3]) +``` + +
+ +## Creating HDF5 files + +ModelArray does not create HDF5 files from raw imaging data — that's the job of the companion conversion tool: + +| Data type | Command | +|:----------|:--------| +| Fixel (`.mif`) | `modelarrayio to-modelarray --backend hdf5 --cohort-file --index-file --directions-file --output ` | +| Voxel (`.nii.gz`) | `modelarrayio to-modelarray --backend hdf5 --cohort-file --mask --output ` | +| Surface (`.dscalar.nii`) | `modelarrayio to-modelarray --backend hdf5 --cohort-file --output ` | + +The command reads the source imaging files listed in a cohort CSV and writes them into the HDF5 layout described above. See the [ModelArrayIO documentation](https://github.com/PennLINC/ModelArrayIO) for detailed usage instructions. + +
+TileDB + +Use the same `modelarrayio to-modelarray` command with `--backend tiledb` and a +`.tdb` output path: + +| Data type | Command | +|:----------|:--------| +| Fixel (`.mif`) | `modelarrayio to-modelarray --backend tiledb --cohort-file --index-file --directions-file --output ` | +| Voxel (`.nii.gz`) | `modelarrayio to-modelarray --backend tiledb --cohort-file --mask --output ` | +| Surface (`.dscalar.nii`) | `modelarrayio to-modelarray --backend tiledb --cohort-file --output ` | + +The command reads the source imaging files listed in a cohort CSV and writes +them into the TileDB-backed layout described above. + +
diff --git a/vignettes/modelling.Rmd b/vignettes/modelling.Rmd index 3a37662..cc3904f 100644 --- a/vignettes/modelling.Rmd +++ b/vignettes/modelling.Rmd @@ -343,24 +343,6 @@ If you want to stream model statistics too (not only transformed scalars), use `write_results_name` and `write_results_file` in `ModelArray.lm()`, `ModelArray.gam()`, or `ModelArray.wrap()`. -#### Quick QA: export one harmonized image to inspect - -After writing `thickness_covfam` into `/scalars`, you can export one row back to an image -and open it in your usual viewer as a spot-check. Use `--column-index` to select which row -from `scalars/thickness_covfam/values` to export. - -```{bash eval=FALSE} -modelarrayio h5-export-nifti-file \ - --input-hdf5 thickness_harmonized.h5 \ - --scalar-name thickness_covfam \ - --column-index 0 \ - --group-mask-file group_mask_thickness.nii.gz \ - --output-file qa_thickness_covfam_col000.nii.gz -``` - -If your workflow uses CIFTI or MIF instead of NIfTI, use the matching export command: -`h5-export-cifti-file` or `h5-export-mif-file`. - ## Modelling across multiple h5 files with `mergeModelArrays()` When scalars live in separate h5 files — for example, cortical thickness in one file @@ -432,13 +414,31 @@ All scalar names must be unique across the inputs. ## Converting results back to image format -After `writeResults()`, use the appropriate ModelArrayIO tool to convert results back to the original image format for visualization: +After `writeResults()`, use the unified `ModelArrayIO` exporter to convert saved +model statistics back to the original image format for visualization. The CLI +entry point is `modelarrayio export-results`, which calls the `export_results()` +exporter and infers the output modality from the arguments you provide. + +For the fixel example used in this vignette: + +```{.console} +$ modelarrayio export-results \ + --index-file ~/Desktop/myProject/FDC/index.mif \ + --directions-file ~/Desktop/myProject/FDC/directions.mif \ + --cohort-file ~/Desktop/myProject/cohort_FDC_n100.csv \ + --analysis-name results_lm \ + --input-hdf5 ~/Desktop/myProject/demo_FDC_n100.h5 \ + --output-dir ~/Desktop/myProject/results_lm +``` + +Use the same command for other ModelArray HDF5 results, changing only the +modality-specific template arguments: -| Data type | Command | Viewer | -|:----------|:--------|:-------| -| Fixel | `fixelstats_write` | MRtrix MRView | -| Voxel | `voxelstats_write` | FSLeyes, MRView | -| Surface | `ciftistats_write` | Connectome Workbench | +| Data type | Command | Extra arguments | Viewer | +|:----------|:--------|:----------------|:-------| +| Fixel | `modelarrayio export-results` | `--index-file`, `--directions-file`, and `--cohort-file` or `--example-file` | MRtrix MRView | +| Voxel | `modelarrayio export-results` | `--mask` | FSLeyes, MRView | +| Surface | `modelarrayio export-results` | `--cohort-file` or `--example-file` | Connectome Workbench | See the [ModelArrayIO documentation](https://github.com/PennLINC/ModelArrayIO) for command-line usage details, and `vignette("walkthrough")` for a worked example with MRView. diff --git a/vignettes/walkthrough.Rmd b/vignettes/walkthrough.Rmd index 2d46a98..751dc6b 100644 --- a/vignettes/walkthrough.Rmd +++ b/vignettes/walkthrough.Rmd @@ -11,6 +11,21 @@ This walkthrough takes you from raw fixel data to statistical results in three s For deeper background on the concepts here, see the [Introductions](../articles/elements.html) vignettes. For more modelling options (GAMs, custom functions), see `vignette("modelling")`. + + ## Step 1. Prepare data and convert to HDF5 ### Download the demo data @@ -39,7 +54,8 @@ This gives you a folder of fixel `.mif` files and a cohort CSV: The CSV file contains one row per subject with columns for covariates and two required columns: - `scalar_name`: the metric name (e.g., `FDC`) -- `source_file`: path to that subject's data file, relative to `--relative-root` +- `source_file`: path to that subject's data file, either absolute or relative + to the directory from which `modelarrayio` is run | subject_id | Age | sex | dti64MeanRelRMS | scalar_name | source_file | |:-----------:|:----:|:---:|:---------------:|:-----------:|:-------------------:| @@ -49,19 +65,49 @@ The CSV file contains one row per subject with columns for covariates and two re ### Convert to HDF5 -Use the `confixel` command from [ModelArrayIO](https://github.com/PennLINC/ModelArrayIO) to convert the `.mif` files into an HDF5 file (see `vignette("hdf5-format")` for what this file contains): +Use the `modelarrayio to-modelarray` command from [ModelArrayIO](https://github.com/PennLINC/ModelArrayIO) to convert the `.mif` files into an HDF5 file (see `vignette("modelarray-storage")` for what this file contains): + +```{.console} +$ conda activate modelarray +$ modelarrayio to-modelarray \ + --backend hdf5 \ + --index-file FDC/index.mif \ + --directions-file FDC/directions.mif \ + --cohort-file cohort_FDC_n100.csv \ + --output demo_FDC_n100.h5 +``` + +
+TileDB + +Use `--backend tiledb` to create a directory-backed TileDB ModelArray store +instead of a single `.h5` file: ```{.console} $ conda activate modelarray -$ confixel \ +$ modelarrayio to-modelarray \ + --backend tiledb \ --index-file FDC/index.mif \ --directions-file FDC/directions.mif \ --cohort-file cohort_FDC_n100.csv \ - --relative-root /home//Desktop/myProject \ - --output-hdf5 demo_FDC_n100.h5 + --output demo_FDC_n100.tdb ``` -For voxel-wise data, use `convoxel` instead. For CIFTI surface data, use `concifti`. See the [ModelArrayIO docs](https://github.com/PennLINC/ModelArrayIO) for details. +TileDB support in R requires the optional `tiledb` and `TileDBArray` packages. +See the optional TileDB installation notes in `vignette("installations")`. + +
+ +Use the same conversion command for other modalities, changing only the +modality-specific arguments: + +| Data type | Command | Extra arguments | +|:----------|:--------|:----------------| +| Fixel | `modelarrayio to-modelarray` | `--index-file` and `--directions-file` | +| Voxel | `modelarrayio to-modelarray` | `--mask` | +| CIFTI | `modelarrayio to-modelarray` | No extra arguments | + +See the [ModelArrayIO docs](https://github.com/PennLINC/ModelArrayIO) for details. ## Step 2. Fit a linear model with ModelArray @@ -83,10 +129,29 @@ modelarray ModelArray located at ~/Desktop/myProject/demo_FDC_n100.h5 Source files: 100 - Scalars: FDC + Scalars: FDC (602229 elements) Analyses: ``` +
+TileDB + +When the path ends in `.tdb`, ModelArray auto-detects the TileDB backend: + +```{r load_tiledb_data, eval=FALSE} +library(ModelArray) + +tdb_path <- "~/Desktop/myProject/demo_FDC_n100.tdb" +csv_path <- "~/Desktop/myProject/cohort_FDC_n100.csv" + +modelarray <- ModelArray(tdb_path, scalar_types = c("FDC")) +modelarray +``` + +You can also be explicit with `backend = "tiledb"` if you prefer. + +
+ ```{r load_csv, eval=FALSE} phenotypes <- read.csv(csv_path) ``` @@ -127,6 +192,24 @@ On a Linux machine with a 10th-gen Xeon at 2.8 GHz using 4 cores, this takes abo writeResults(h5_path, df.output = result_full, analysis_name = "results_lm") ``` +
+TileDB + +Write the same result table back into a TileDB store by using the `.tdb` path. +The explicit backend argument is optional for `.tdb` paths, but shown here for +clarity: + +```{r write_tiledb_results, eval=FALSE} +writeResults( + fn.output = tdb_path, + df.output = result_full, + analysis_name = "results_lm", + backend = "tiledb" +) +``` + +
+ You can verify the results were saved: ```{r verify_results, eval=FALSE} @@ -141,29 +224,81 @@ modelarray_new ModelArray located at ~/Desktop/myProject/demo_FDC_n100.h5 Source files: 100 - Scalars: FDC + Scalars: FDC (602229 elements) Analyses: results_lm ``` +
+TileDB + +Load saved analyses from the TileDB store the same way: + +```{r verify_tiledb_results, eval=FALSE} +modelarray_new <- ModelArray( + tdb_path, + scalar_types = "FDC", + analysis_names = "results_lm" +) +modelarray_new +``` + +
+ ## Step 3. Convert results back and visualize ### Convert to fixel `.mif` format -Use the `fixelstats_write` command from `ModelArrayIO` to convert the results back to `.mif` files for visualization: +Use the `modelarrayio export-results` command from `ModelArrayIO` to convert the results back to `.mif` files for visualization: ```{.console} $ conda activate modelarray -$ fixelstats_write \ +$ modelarrayio export-results \ --index-file FDC/index.mif \ --directions-file FDC/directions.mif \ --cohort-file cohort_FDC_n100.csv \ - --relative-root /home//Desktop/myProject \ --analysis-name results_lm \ --input-hdf5 demo_FDC_n100.h5 \ --output-dir results_lm ``` -For voxel data, use `voxelstats_write`. For CIFTI, use `ciftistats_write`. +
+TileDB + +`modelarrayio export-results` currently only exports statistical results from HDF5 +ModelArray files via `--input-hdf5`. If you model with TileDB and need viewer +files, keep or write an HDF5 results file for the export step: + +```{r write_hdf5_export_copy, eval=FALSE} +writeResults( + fn.output = "demo_FDC_n100_results_for_export.h5", + df.output = result_full, + analysis_name = "results_lm", + backend = "hdf5" +) +``` + +Then export from that HDF5 results file: + +```{.console} +$ modelarrayio export-results \ + --index-file FDC/index.mif \ + --directions-file FDC/directions.mif \ + --cohort-file cohort_FDC_n100.csv \ + --analysis-name results_lm \ + --input-hdf5 demo_FDC_n100_results_for_export.h5 \ + --output-dir results_lm +``` + +
+ +Use the same export command for other ModelArray HDF5 results, changing only the +modality-specific template arguments: + +| Data type | Command | Extra arguments | Viewer | +|:----------|:--------|:----------------|:-------| +| Fixel | `modelarrayio export-results` | `--index-file`, `--directions-file`, and `--cohort-file` or `--example-file` | MRtrix MRView | +| Voxel | `modelarrayio export-results` | `--mask` | FSLeyes, MRView | +| CIFTI | `modelarrayio export-results` | `--cohort-file` or `--example-file` | Connectome Workbench | ### View in MRView @@ -184,6 +319,6 @@ $ mrview ## Next steps - **More model types**: fit GAMs with nonlinear smooth terms, or run arbitrary custom functions with `ModelArray.wrap()` — see `vignette("modelling")` -- **Explore your data**: use ModelArray's accessor functions to inspect scalars, results, and per-element data frames — see `vignette("exploring-h5")` +- **Explore your data**: use ModelArray's accessor functions to inspect scalars, results, and per-element data frames — see `vignette("exploring-modelarray-data")` - **Scale up**: split large analyses across HPC jobs — see `vignette("element-splitting")` -- **Understand the file format**: learn about chunking, compression, and memory efficiency — see `vignette("hdf5-large-analyses")` +- **Understand large-scale analyses**: learn about chunking, compression, and memory efficiency — see `vignette("large-scale-analyses")`