diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6bb84762f..31021044f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,7 +28,10 @@ jobs: run: | git config --global user.email "runner@github.com" git config --global user.name "GitHub Runner" - - uses: actions/checkout@v3 + - name: Checkout repository + uses: actions/checkout@v3 + with: + submodules: true - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v4 with: diff --git a/.github/workflows/docs-preview.yml b/.github/workflows/docs-preview.yml index 9025eb8af..8c98b4c35 100644 --- a/.github/workflows/docs-preview.yml +++ b/.github/workflows/docs-preview.yml @@ -18,12 +18,11 @@ concurrency: preview-${{ github.ref }} jobs: build-docs: runs-on: ubuntu-latest - strategy: - fail-fast: false - steps: - name: Checkout source uses: actions/checkout@v3 + with: + submodules: true - name: Set up Python 3.10 uses: actions/setup-python@v4 with: @@ -51,4 +50,4 @@ jobs: if: github.event_name == 'pull_request' uses: rossjrw/pr-preview-action@v1 with: - source-dir: docs/_build/main \ No newline at end of file + source-dir: docs/_build/main diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index e7ed4c66a..2b27452c1 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -12,15 +12,13 @@ on: jobs: build-docs: runs-on: ubuntu-latest - strategy: - fail-fast: false - steps: - name: Checkout source uses: actions/checkout@v3 with: # require all of history to see all tagged versions' docs fetch-depth: 0 + submodules: true - name: Set up Python 3.10 uses: actions/setup-python@v4 with: diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 000000000..0c541549d --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "junifer/external/h5io"] + path = junifer/external/h5io + url = https://github.com/juaml/h5io diff --git a/docs/changes/latest.inc b/docs/changes/latest.inc index ff0da3932..23c012ce6 100644 --- a/docs/changes/latest.inc +++ b/docs/changes/latest.inc @@ -46,6 +46,9 @@ Enhancements - Add :class:`junifer.markers.TemporalSNRParcels` and :class:`junifer.markers.TemporalSNRSpheres` (:gh:`163` by `Leonard Sasse`_). +- Add support for HDF5 feature storage via :class:`junifer.storage.HDF5FeatureStorage` + (:gh:`147` by `Synchon Mandal`_). + Bugs ~~~~ diff --git a/docs/conf.py b/docs/conf.py index a11f4ed53..a72c37e38 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -26,7 +26,7 @@ sys.path.append((curdir / "sphinxext").as_posix()) # -- Project information ----------------------------------------------------- project = "junifer" -copyright = "2022, Authors of junifer" +copyright = "2023, Authors of junifer" author = "Fede Raimondo" # -- General configuration --------------------------------------------------- diff --git a/docs/understanding/storage.rst b/docs/understanding/storage.rst index 5c8d68495..fb637c91b 100644 --- a/docs/understanding/storage.rst +++ b/docs/understanding/storage.rst @@ -67,4 +67,8 @@ Currently supported storage interfaces * - :class:`junifer.storage.SQLiteFeatureStorage` - ``.sqlite`` - SQLite - - ``matrix``, ``table``, ``timeseries`` + - ``matrix``, ``vector``, ``timeseries`` + * - :class:`junifer.storage.HDF5FeatureStorage` + - ``.hdf5`` + - HDF5 + - ``matrix``, ``vector``, ``timeseries`` diff --git a/junifer/datagrabber/aomic/id1000.py b/junifer/datagrabber/aomic/id1000.py index 1a1b188e2..428192b46 100644 --- a/junifer/datagrabber/aomic/id1000.py +++ b/junifer/datagrabber/aomic/id1000.py @@ -7,7 +7,7 @@ # License: AGPL from pathlib import Path -from typing import Union, Dict +from typing import Dict, Union from junifer.datagrabber import PatternDataladDataGrabber diff --git a/junifer/datagrabber/aomic/piop2.py b/junifer/datagrabber/aomic/piop2.py index b9fad565e..284d259a3 100644 --- a/junifer/datagrabber/aomic/piop2.py +++ b/junifer/datagrabber/aomic/piop2.py @@ -7,7 +7,7 @@ # License: AGPL from pathlib import Path -from typing import List, Union, Dict +from typing import Dict, List, Union from junifer.datagrabber import PatternDataladDataGrabber diff --git a/junifer/external/h5io b/junifer/external/h5io new file mode 160000 index 000000000..1d2ba9925 --- /dev/null +++ b/junifer/external/h5io @@ -0,0 +1 @@ +Subproject commit 1d2ba9925c12f5e3c191eca888b829e43c206658 diff --git a/junifer/storage/__init__.py b/junifer/storage/__init__.py index b8c4f8942..1e756e6cc 100644 --- a/junifer/storage/__init__.py +++ b/junifer/storage/__init__.py @@ -7,3 +7,4 @@ from .base import BaseFeatureStorage from .pandas_base import PandasBaseFeatureStorage from .sqlite import SQLiteFeatureStorage +from .hdf5 import HDF5FeatureStorage diff --git a/junifer/storage/base.py b/junifer/storage/base.py index 5d2aba2e4..04cdcb71f 100644 --- a/junifer/storage/base.py +++ b/junifer/storage/base.py @@ -190,7 +190,7 @@ class BaseFeatureStorage(ABC): data: np.ndarray, col_names: Optional[Iterable[str]] = None, row_names: Optional[Iterable[str]] = None, - matrix_kind: Optional[str] = "full", + matrix_kind: str = "full", diagonal: bool = True, ) -> None: """Store matrix. diff --git a/junifer/storage/hdf5.py b/junifer/storage/hdf5.py new file mode 100644 index 000000000..7ccc1fcac --- /dev/null +++ b/junifer/storage/hdf5.py @@ -0,0 +1,921 @@ +"""Provide concrete implementation for feature storage via HDF5.""" + +# Authors: Synchon Mandal +# Federico Raimondo +# License: AGPL + + +from collections import defaultdict +from functools import reduce +from pathlib import Path +from typing import Any, Dict, Iterable, List, Optional, Union + +import numpy as np +import pandas as pd +from tqdm import tqdm + +from ..api.decorators import register_storage +from ..external.h5io.h5io import ChunkedArray, read_hdf5, write_hdf5 +from ..utils import logger, raise_error +from .base import BaseFeatureStorage +from .utils import element_to_prefix, matrix_to_vector, store_matrix_checks + + +@register_storage +class HDF5FeatureStorage(BaseFeatureStorage): + """Concrete implementation for feature storage via HDF5. + + Parameters + ---------- + uri : str or pathlib.Path + The path to the file to be used. + single_output : bool, optional + If False, will create one HDF5 file per element. The name + of the file will be prefixed with the respective element. + If True, will create only one HDF5 file as specified in the + ``uri`` and store all the elements in the same file. Concurrent + writes should be handled with care (default True). + overwrite : bool or "update", optional + Whether to overwrite existing file. If True, will overwrite and + if "update", will update existing entry or append (default "update"). + compression : {0-9}, optional + Level of gzip compression: 0 (lowest) to 9 (highest) (default 7). + force_float32 : bool, optional + Whether to force casting of numpy.ndarray values to float32 if float64 + values are found (default True). + chunk_size : int, optional + The chunk size to use when collecting data from element files in + :meth:`junifer.storage.HDF5FeatureStorage.collect`. If the file count + is smaller than the value, the minimum is used (default 100). + + See Also + -------- + SQLiteFeatureStorage : The concrete class for SQLite-based feature storage. + + """ + + def __init__( + self, + uri: Union[str, Path], + single_output: bool = True, + overwrite: Union[bool, str] = "update", + compression: int = 7, + force_float32: bool = True, + chunk_size: int = 100, + ) -> None: + # Convert str to Path + if not isinstance(uri, Path): + uri = Path(uri) + + # Create parent directories if not present + if not uri.parent.exists(): + logger.info( + f"Output directory: '{uri.parent.resolve()}' " + "does not exist, creating now" + ) + uri.parent.mkdir(parents=True, exist_ok=True) + + # Available storage kinds + storage_types = ["vector", "timeseries", "matrix"] + + super().__init__( + uri=uri, + storage_types=storage_types, + single_output=single_output, + ) + + self.overwrite = overwrite + self.compression = compression + self.force_float32 = force_float32 + self.chunk_size = chunk_size + + def get_valid_inputs(self) -> List[str]: + """Get valid storage types for input. + + Returns + ------- + list of str + The list of storage types that can be used as input for this + storage. + + """ + return ["matrix", "vector", "timeseries"] + + def _fetch_correct_uri_for_io(self, element: Optional[Dict]) -> str: + """Return proper URI for I/O based on `element`. + + If `element` is None, will return `self.uri`. + + Parameters + ---------- + element : dict, optional + The element as dictionary (default None). + + Returns + ------- + str + Formatted URI for accessing metadata and data. + + """ + if not self.single_output and not element: + raise_error( + msg=( + "`element` must be provided when `single_output` " + "is False" + ), + klass=RuntimeError, + ) + elif not self.single_output and element: + # element access for multi output only + prefix = element_to_prefix(element=element) + else: + # parent access for single output, ignore element + prefix = "" + # Format URI based on prefix + return f"{self.uri.parent}/{prefix}{self.uri.name}" # type: ignore + + def _read_metadata( + self, element: Optional[Dict[str, str]] = None + ) -> Dict[str, Dict[str, Any]]: + """Read metadata (should not be called directly). + + Parameters + ---------- + element : dict, optional + The element as dictionary (default None). + + Returns + ------- + dict of dict + The stored metadata for the element. + + Raises + ------ + IOError + If HDF5 file or `meta` does not exist. + + """ + # Get correct URI for element; + # is different from uri if single_output is False + uri = self._fetch_correct_uri_for_io(element=element) + + try: + logger.info(f"Loading HDF5 metadata from: {uri}") + metadata = read_hdf5( + fname=uri, + title="meta", + slash="ignore", + ) + except IOError: + raise_error( + msg=f"HDF5 file not found at: {uri}", + klass=IOError, + ) + except ValueError: + raise_error( + msg=f"`meta` not found in: {uri}", + klass=IOError, + ) + else: + logger.info(f"Loaded HDF5 metadata from: {uri}") + return metadata + + def list_features(self) -> Dict[str, Dict[str, Any]]: + """List the features in the storage. + + Returns + ------- + dict + List of features in the storage. The keys are the feature MD5 to + be used in :meth:`junifer.storage.HDF5FeatureStorage.read_df` + and the values are the metadata of each feature. + + """ + # Read metadata + metadata = read_hdf5( + fname=str(self.uri.resolve()), # type: ignore + title="meta", + slash="ignore", + ) + return metadata + + def _read_data( + self, md5: str, element: Optional[Dict[str, str]] = None + ) -> Dict[str, Any]: + """Read data (should not be called directly). + + Parameters + ---------- + md5 : str + The MD5 used as the HDF5 group name. + element : dict, optional + The element as dictionary (default None). + + Returns + ------- + dict + The retrieved data. + + Raises + ------ + IOError + If HDF5 file or data does not exist. + + """ + # Get correct URI for element; + # is different from uri if single_output is False + uri = self._fetch_correct_uri_for_io(element=element) + + try: + logger.info(f"Loading HDF5 data for {md5} from: {uri}") + data = read_hdf5( + fname=uri, + title=md5, + slash="ignore", + ) + except IOError: + raise_error( + msg=f"HDF5 file not found at: {uri}", + klass=IOError, + ) + except ValueError: + raise_error( + msg=f"`{md5}` not found in: {uri}", + klass=IOError, + ) + else: + logger.info(f"Loaded HDF5 data for {md5} from: {uri}") + return data + + def read_df( + self, + feature_name: Optional[str] = None, + feature_md5: Optional[bool] = None, + ) -> pd.DataFrame: + """Read feature into a pandas.DataFrame. + + Either one of ``feature_name`` or ``feature_md5`` needs to be + specified. + + Parameters + ---------- + feature_name : str, optional + Name of the feature to read (default None). + feature_md5 : str, optional + MD5 hash of the feature to read (default None). + + Returns + ------- + pandas.DataFrame + The features as a dataframe. + + Raises + ------ + IOError + If HDF5 file does not exist. + + """ + # Parameter conflict + if feature_md5 and feature_name: + raise_error( + msg=( + "Only one of `feature_name` or `feature_md5` can be " + "specified." + ) + ) + # Parameter absence + elif not feature_md5 and not feature_name: + raise_error( + msg=( + "At least one of `feature_name` or `feature_md5` " + "must be specified." + ) + ) + # Parameter check pass; read metadata + metadata = read_hdf5( + fname=str(self.uri.resolve()), # type: ignore + title="meta", + slash="ignore", + ) + # Initialize MD5 variable + md5: str = "" + + # Consider feature_md5 + if feature_md5: + logger.debug( + f"Validating feature MD5 '{feature_md5}' in metadata " + f"for: {self.uri.resolve()} ..." # type: ignore + ) + # Validate MD5 + if feature_md5 in metadata: + md5 = feature_md5 # type: ignore + else: + raise_error(msg=f"Feature MD5 '{feature_md5}' not found") + + # Consider feature_name + elif feature_name: + logger.debug( + f"Validating feature name '{feature_name}' in metadata " + f"for: {self.uri.resolve()} ..." # type: ignore + ) + # Retrieve MD5 for feature_name + # Implicit counter for duplicate feature_name with different + # MD5; happens when something is wrong with marker computation + feature_name_duplicates_with_different_md5 = [] + for md5, meta in metadata.items(): + if meta["name"] == feature_name: + feature_name_duplicates_with_different_md5.append(md5) + + # Check for no / duplicate feature_name + if len(feature_name_duplicates_with_different_md5) == 0: + raise_error(msg=f"Feature '{feature_name}' not found") + elif len(feature_name_duplicates_with_different_md5) > 1: + raise_error( + msg=( + "More than one feature with name " + f"'{feature_name}' found. You can bypass this " + "issue by specifying a `feature_md5`." + ) + ) + + md5 = feature_name_duplicates_with_different_md5[0] + + # Read data from HDF5 + hdf_data = read_hdf5( + fname=str(self.uri.resolve()), # type: ignore + title=md5, + slash="ignore", + ) + reshaped_data = hdf_data["data"] + + # Generate index for the data + logger.debug(f"Generating pandas.MultiIndex for {md5} ...") + + if hdf_data["kind"] == "matrix": + # Set index for element + element_idx = hdf_data["element"] + # Flatten data and get column headers for dataframe + flat_data, columns = matrix_to_vector( + data=hdf_data["data"], + col_names=hdf_data["column_headers"], + row_names=hdf_data["row_headers"], + matrix_kind=hdf_data["matrix_kind"], + diagonal=bool(hdf_data["diagonal"]), + ) + # Convert data to proper 2D + reshaped_data = flat_data.T + elif hdf_data["kind"] == "vector": + # Set index for element + element_idx = hdf_data["element"] + # Set column headers for dataframe + columns = hdf_data["column_headers"] + # Convert data to proper 2D + reshaped_data = hdf_data["data"].T + elif hdf_data["kind"] == "timeseries": + # Create dictionary for aggregating index data + element_idx = defaultdict(list) + for idx, element in enumerate(hdf_data["element"]): + # Get row count for the element + n_rows, _ = hdf_data["data"][:, :, idx].shape + # Set rows for the index + for key, val in element.items(): + element_idx[key].extend([val] * n_rows) + # Add extra column for timepoints + element_idx[hdf_data["row_header_column_name"]].extend( + np.arange(n_rows) + ) + # Set column headers for dataframe + columns = hdf_data["column_headers"] + # Convert data from 3D to 2D + reshaped_data = hdf_data["data"].reshape(-1, 1) + + # Create dataframe for index + idx_df = pd.DataFrame(data=element_idx) # type: ignore + # Create multiindex from dataframe + hdf_data_idx = pd.MultiIndex.from_frame(df=idx_df) + logger.debug(f"Generated pandas.MultiIndex for {md5} ...") + + # Convert to DataFrame + logger.debug(f"Converting HDF5 data for {md5} to pandas.DataFrame ...") + df = pd.DataFrame( + data=reshaped_data, + index=hdf_data_idx, + columns=columns, # type: ignore + dtype=hdf_data["data"].dtype, + ) + logger.debug(f"Converted HDF5 data for {md5} to pandas.DataFrame ...") + + return df + + def _write_processed_data( + self, fname: str, processed_data: Dict[str, Any], title: str + ) -> None: + """Write processed data to HDF5 (should not be called directly). + + This is used primarily in + :func:`junifer.storage.HDF5FeatureStorage.store_metadata` and + ``_store_data``. + + Parameters + ---------- + fname : str + The HDF5 file to store data. + processed_data : dict + The processed data as dictionary. + title : str + The top-level directory name. + + """ + logger.debug(f"Writing processed HDF5 data to: {fname} ...") + # Write to HDF5 + write_hdf5( + fname=fname, + data=processed_data, + overwrite=self.overwrite, # type: ignore + compression=self.compression, + title=title, + slash="error", + use_json=True, + ) + logger.debug(f"Wrote processed HDF5 data to: {fname} ...") + + def store_metadata( + self, + meta_md5: str, + element: Dict[str, str], + meta: Dict[str, Any], + ) -> None: + """Store metadata. + + This method first loads existing metadata (if any) using + ``_read_metadata`` and appends to it the new metadata and then saves + the updated metadata using ``_write_processed_data``. It will only + store metadata if ``meta_md5`` is not found already. + + Parameters + ---------- + meta_md5 : str + The metadata MD5 hash. + element : dict + The element as a dictionary. + meta : dict + The metadata as a dictionary. + + """ + # Read metadata; if no file found, create an empty dictionary + try: + metadata = self._read_metadata(element=element) + except IOError: + logger.debug(f"Creating new metadata map for {meta_md5} ...") + metadata = {} + + # Only add entry if MD5 is not present + if meta_md5 not in metadata: + logger.debug(f"HDF5 metadata for {meta_md5} not found, adding ...") + # Update metadata + metadata[meta_md5] = meta + + # Get correct URI for element; + # is different from uri if single_output is False + uri = self._fetch_correct_uri_for_io(element=element) + + logger.info(f"Writing HDF5 metadata for {meta_md5} to: {uri}") + logger.debug(f"HDF5 overwrite is set to: {self.overwrite} ...") + logger.debug( + "HDF5 gzip compression level is set to: " + f"{self.compression} ..." + ) + + # Write metadata + self._write_processed_data( + fname=uri, + processed_data=metadata, + title="meta", + ) + + logger.info(f"Wrote HDF5 metadata for {meta_md5} to: {uri}") + else: + logger.debug( + f"HDF5 metadata for {meta_md5} found, skipping store ..." + ) + + def _store_data( + self, + kind: str, + meta_md5: str, + element: List[Dict[str, str]], + data: np.ndarray, + **kwargs: Any, + ) -> None: + """Store data. + + This method first loads existing data (if any) using + ``_read_data`` and appends to it the `element` and `data` + values, and writes them and other information passed via + `**kwargs` using ``_write_processed_data``. + + Parameters + ---------- + kind : {"matrix", "vector", "timeseries"} + The storage kind. + meta_md5 : str + The metadata MD5 hash. + element : list of dict + The element as list of dictionary. + data : numpy.ndarray + The data to store. + **kwargs : dict + Keyword arguments passed from the calling method. + + """ + # Read existing data; if no file found, create an empty dictionary + try: + stored_data = self._read_data(md5=meta_md5, element=element[0]) + except IOError: + logger.debug(f"Creating new data map for {meta_md5} ...") + stored_data = {} + + # Initialize dictionary to aggregate data to write + data_to_write = kwargs + + # Optional casting of float64 values to float32 for numpy.ndarray + if data.dtype == np.dtype("float64") and self.force_float32: + data = data.astype(dtype=np.dtype("float32"), casting="same_kind") + + # Handle cases for existing and new entry + if not stored_data: + logger.debug(f"Writing new data for {meta_md5} ...") + # New entry; add as is + data_to_write.update( + { + "element": element, + "data": data, + # for serialization / deserialization of storage type + "kind": kind, + } + ) + elif stored_data: + # Set up stored kwargs + stored_kwargs = [ + key + for key in stored_data.keys() + if key not in ("element", "data") + ] + # Set up to be stored kwargs + to_be_stored_kwargs = kwargs + # Update with kind + to_be_stored_kwargs["kind"] = kind + # Validate the kwargs + if set(stored_kwargs) != set(to_be_stored_kwargs): + raise_error( + msg=( + f"The additional data for {meta_md5} do not match " + "the ones already stored. This can be due " + "to some changes in marker computation, please " + "verify and try again." + ), + klass=RuntimeError, + ) + + # Check for duplicate elements; if found, return immediately + logger.debug(f"Checking duplicate elements for {meta_md5} ...") + for stored_element in stored_data["element"]: + if stored_element == element[0]: + logger.info( + f"Duplicate element: {element[0]} found for " + f"{meta_md5}, skipping store ... " + ) + return None + logger.debug(f"No duplicate elements found for {meta_md5} ...") + + logger.debug( + f"Existing data found for {meta_md5}, appending to it ..." + ) + # Existing entry; append to existing + # "element" and "data" + data_to_write.update( + { + "element": stored_data["element"] + element, + "data": np.concatenate( + (stored_data["data"], data), axis=-1 + ), + # for serialization / deserialization of storage type + "kind": kind, + } + ) + + # Get correct URI for element; + # is different from uri if single_output is False + uri = self._fetch_correct_uri_for_io(element=element[0]) + + logger.info(f"Writing HDF5 data for {meta_md5} to: {uri}") + logger.debug(f"HDF5 overwrite is set to: {self.overwrite} ...") + logger.debug( + f"HDF5 gzip compression level is set to: {self.compression} ..." + ) + + # Write data + self._write_processed_data( + fname=uri, + processed_data=data_to_write, + title=meta_md5, + ) + + logger.info(f"Wrote HDF5 data for {meta_md5} to: {uri}") + + def store_matrix( + self, + meta_md5: str, + element: Dict[str, str], + data: np.ndarray, + col_names: Optional[Iterable[str]] = None, + row_names: Optional[Iterable[str]] = None, + matrix_kind: str = "full", + diagonal: bool = True, + row_header_col_name: str = "ROI", + ) -> None: + """Store matrix. + + This method performs parameter checks and then calls + ``_store_data`` for storing the data. + + Parameters + ---------- + meta_md5 : str + The metadata MD5 hash. + element : dict + The element as dictionary. + data : numpy.ndarray + The matrix data to store. + col_names : list or tuple of str, optional + The column labels (default None). + row_names : list or tuple of str, optional + The row labels (default None). + matrix_kind : str, optional + The kind of matrix: + + * ``triu`` : store upper triangular only + * ``tril`` : store lower triangular + * ``full`` : full matrix + + (default "full"). + diagonal : bool, optional + Whether to store the diagonal. If ``matrix_kind`` is "full", + setting this to False will raise an error (default True). + row_header_col_name : str, optional + The column name for the row header column (default "ROI"). + + Raises + ------ + ValueError + If invalid ``matrix_kind`` is provided, ``diagonal = False`` + for ``matrix_kind = "full"``, non-square data is provided + for ``matrix_kind = {"triu", "tril"}``, length of ``row_names`` + do not match data row count, or length of ``col_names`` do not + match data column count. + + """ + # Row data validation + if row_names is None: + row_names = [f"r{i}" for i in range(data.shape[0])] + # Column data validation + if col_names is None: + col_names = [f"c{i}" for i in range(data.shape[1])] + # Parameter checks + store_matrix_checks( + matrix_kind=matrix_kind, + diagonal=diagonal, + data_shape=data.shape, + row_names_len=len(row_names), # type: ignore + col_names_len=len(col_names), # type: ignore + ) + # Store + self._store_data( + kind="matrix", + meta_md5=meta_md5, + element=[element], # convert to list + data=data[:, :, np.newaxis], # convert to 3D + column_headers=col_names, + row_headers=row_names, + matrix_kind=matrix_kind, + diagonal=diagonal, + row_header_column_name=row_header_col_name, + ) + + def store_vector( + self, + meta_md5: str, + element: Dict[str, str], + data: Union[np.ndarray, List], + col_names: Optional[Iterable[str]] = None, + ) -> None: + """Store vector. + + Parameters + ---------- + meta_md5 : str + The metadata MD5 hash. + element : dict + The element as dictionary. + data : numpy.ndarray or list + The vector data to store. + col_names : list or tuple of str, optional + The column labels (default None). + + """ + if isinstance(data, list): + logger.debug( + f"Flattening and converting vector data list for {meta_md5}, " + "to numpy.ndarray ..." + ) + # Flatten out list and convert to np.ndarray + processed_data = np.array(np.ravel(data)) + elif isinstance(data, np.ndarray): + logger.debug( + f"Flattening vector data numpy.ndarray for {meta_md5} ..." + ) + # Flatten out array + processed_data = data.ravel() + + self._store_data( + kind="vector", + meta_md5=meta_md5, + element=[element], # convert to list + data=processed_data[:, np.newaxis], # convert to 2D + column_headers=col_names, + ) + + def store_timeseries( + self, + meta_md5: str, + element: Dict[str, str], + data: np.ndarray, + col_names: Optional[Iterable[str]] = None, + ) -> None: + """Store timeseries. + + Parameters + ---------- + meta_md5 : str + The metadata MD5 hash. + element : dict + The element as dictionary. + data : numpy.ndarray + The timeseries data to store. + col_names : list or tuple of str, optional + The column labels (default None). + + """ + self._store_data( + kind="timeseries", + meta_md5=meta_md5, + element=[element], # convert to list + data=data[:, :, np.newaxis], # convert to 3D + column_headers=col_names, + row_header_column_name="timepoint", + ) + + def collect(self) -> None: + """Implement data collection. + + This method globs the element files and runs a loop + over them while reading metadata and then runs a loop + over all the stored features in the metadata, storing + the metadata and the feature data right after reading. + + Raises + ------ + NotImplementedError + If ``single_output`` is True. + + """ + if self.single_output is True: + raise_error( + msg="collect() is not implemented for single output.", + klass=NotImplementedError, + ) + + # Glob files + globbed_files = self.uri.parent.glob( # type: ignore + f"*{self.uri.name}" # type: ignore + ) + + # Create new storage instance + out_storage = HDF5FeatureStorage(uri=self.uri, overwrite="update") + + # Run loop to collect metadata + logger.info( + "Collecting metadata from " + f"{self.uri.parent}/*{self.uri.name}" # type: ignore + ) + # Collect element files per feature MD5 + elements_per_feature_md5 = defaultdict(list) + for file_ in tqdm(globbed_files, desc="file-metadata"): + logger.debug(f"Reading HDF5 file: {file_} ...") + # Create new storage instance to load metadata + in_storage = HDF5FeatureStorage(uri=file_) + # Load metadata from new instance + in_metadata = in_storage._read_metadata() + + logger.info(f"Updating HDF5 metadata with metadata from: {file_}") + # Load metadata; empty dictionary if first entry; + # can be replaced with store_metadata() if run on a loop + # for the metadata entries from in_storage + try: + out_metadata = out_storage._read_metadata() + except IOError: + out_metadata = {} + # Update metadata + out_metadata.update(in_metadata) + # Save metadata + out_storage._write_processed_data( + fname=str(self.uri.resolve()), # type: ignore + processed_data=out_metadata, + title="meta", + ) + # Update element files for found MD5s + for feature_md5 in in_metadata.keys(): + elements_per_feature_md5[feature_md5].append(file_) + + # Run loop to collect data per feature per file + logger.info( + "Collecting data from " + f"{self.uri.parent}/*{self.uri.name}" # type: ignore + ) + for feature_md5, element_files in tqdm( + elements_per_feature_md5.items(), desc="feature" + ): + element_count = len(element_files) + # Chunk size for collecting + chunk_size = min(self.chunk_size, element_count) + # Operate on chunks + for chunk_idx, chunk_start in tqdm( + enumerate(range(0, element_count, chunk_size)), desc="chunk" + ): + # Store the chunk files' data + stored_data_for_chunk: List[Dict[str, Any]] = [] + # Read the files of a chunk + for i in tqdm( + range(chunk_start, chunk_start + chunk_size), + desc="file-data", + ): + file_ = element_files[i] + logger.debug( + f"Reading feature MD5: '{feature_md5}' " + f"from HDF5 file: {file_} ..." + ) + # Read from HDF5 and collect data + stored_data_for_chunk.append( + read_hdf5( + fname=str(file_), + title=feature_md5, + slash="ignore", + ) + ) + + # Concatenate the features data for a chunk + features_data = np.concatenate( + [x["data"] for x in stored_data_for_chunk], axis=-1 + ) + # Make dictionary to write the collected data; + # first the static data then the dynamic data + data_to_write = { + key: val + for key, val in stored_data_for_chunk[0].items() + if key not in ("data", "element") + } + # Join the features element for a chunk + data_to_write["element"] = reduce( + lambda acc, elem: acc + elem, + [x["element"] for x in stored_data_for_chunk], + [], + ) + # Write data in chunks to avoid memory usage spikes + # Start with the case for 2D + array_shape = [features_data.shape[0]] + array_chunk_size = [features_data.shape[0]] + # Append second dimension for 3D + if features_data.ndim == 3: + array_shape.append(features_data.shape[1]) + array_chunk_size.append(features_data.shape[1]) + # Append final dimension of element count + array_shape.append(element_count) + # Append final dimension of chunk size + array_chunk_size.append(chunk_size) + # Write chunked array + data_to_write["data"] = ChunkedArray( + data=features_data, + shape=tuple(array_shape), + chunk_size=tuple(array_chunk_size), + n_chunk=chunk_idx, + ) + # Write to HDF5 + write_hdf5( + fname=str(self.uri.resolve()), # type: ignore + data=data_to_write, + overwrite=self.overwrite, # type: ignore + compression=0, + title=feature_md5, + slash="error", + use_json=False, + ) diff --git a/junifer/storage/sqlite.py b/junifer/storage/sqlite.py index 706dc64ae..e511c9e49 100644 --- a/junifer/storage/sqlite.py +++ b/junifer/storage/sqlite.py @@ -18,7 +18,7 @@ from tqdm import tqdm from ..api.decorators import register_storage from ..utils import logger, raise_error, warn_with_log from .pandas_base import PandasBaseFeatureStorage -from .utils import element_to_prefix +from .utils import element_to_prefix, matrix_to_vector, store_matrix_checks if TYPE_CHECKING: @@ -49,6 +49,7 @@ class SQLiteFeatureStorage(PandasBaseFeatureStorage): See Also -------- PandasBaseFeatureStorage : The base class for Pandas-based feature storage. + HDF5FeatureStorage : The concrete class for HDF5-based feature storage. """ @@ -399,7 +400,7 @@ class SQLiteFeatureStorage(PandasBaseFeatureStorage): data: np.ndarray, col_names: Optional[List[str]] = None, row_names: Optional[List[str]] = None, - matrix_kind: Optional[str] = "full", + matrix_kind: str = "full", diagonal: bool = True, ) -> None: """Store matrix. @@ -431,56 +432,28 @@ class SQLiteFeatureStorage(PandasBaseFeatureStorage): this to False will raise an error (default True). """ - if diagonal is False and matrix_kind not in ["triu", "tril"]: - raise_error( - msg="Diagonal cannot be False if kind is not full", - klass=ValueError, - ) - - if matrix_kind in ["triu", "tril"]: - if data.shape[0] != data.shape[1]: - raise_error( - "Cannot store a non-square matrix as a triangular matrix", - klass=ValueError, - ) - - if matrix_kind == "triu": - k = 0 if diagonal is True else 1 - data_idx = np.triu_indices(data.shape[0], k=k) - elif matrix_kind == "tril": - k = 0 if diagonal is True else -1 - data_idx = np.tril_indices(data.shape[0], k=k) - elif matrix_kind == "full": - data_idx = ( - np.repeat(np.arange(data.shape[0]), data.shape[1]), - np.tile(np.arange(data.shape[1]), data.shape[0]), - ) - else: - raise_error(msg=f"Invalid kind {matrix_kind}", klass=ValueError) - + # Row data validation if row_names is None: row_names = [f"r{i}" for i in range(data.shape[0])] - elif len(row_names) != data.shape[0]: - raise_error( - msg="Number of row names does not match number of rows", - klass=ValueError, - ) - + # Column data validation if col_names is None: col_names = [f"c{i}" for i in range(data.shape[1])] - elif len(col_names) != data.shape[1]: - raise_error( - msg="Number of column names does not match number of columns", - klass=ValueError, - ) - - # Subset data - flat_data = data[data_idx] - # Generate flat 1D row X column names - columns = [ - f"{row_names[i]}~{col_names[j]}" - for i, j in zip(data_idx[0], data_idx[1]) - ] + # Parameter checks + store_matrix_checks( + matrix_kind=matrix_kind, + diagonal=diagonal, + data_shape=data.shape, + row_names_len=len(row_names), # type: ignore + col_names_len=len(col_names), # type: ignore + ) + # Matrix to vector conversion + flat_data, columns = matrix_to_vector( + data=data, + col_names=col_names, + row_names=row_names, + matrix_kind=matrix_kind, + diagonal=diagonal, + ) # Convert element metadata to index n_rows = 1 @@ -488,7 +461,9 @@ class SQLiteFeatureStorage(PandasBaseFeatureStorage): element=element, n_rows=n_rows, rows_col_name=None ) # Prepare new dataframe - data_df = pd.DataFrame(flat_data[None, :], columns=columns, index=idx) + data_df = pd.DataFrame( + flat_data[np.newaxis, :], columns=columns, index=idx + ) # SQLite's SQLITE_MAX_COLUMN is 2000, so if more than that, # convert it to long format diff --git a/junifer/storage/tests/test_hdf5.py b/junifer/storage/tests/test_hdf5.py new file mode 100644 index 000000000..a78fb705d --- /dev/null +++ b/junifer/storage/tests/test_hdf5.py @@ -0,0 +1,948 @@ +"""Provide tests for HDF5 storage interface.""" + +# Authors: Synchon Mandal +# Federico Raimondo +# License: AGPL + +from pathlib import Path + +import numpy as np +import pytest +from numpy.testing import assert_array_equal +from pandas.testing import assert_frame_equal + +from junifer.storage import HDF5FeatureStorage +from junifer.storage.utils import element_to_prefix, process_meta + + +def test_get_valid_inputs() -> None: + """Test valid inputs.""" + storage = HDF5FeatureStorage(uri="/tmp") + assert storage.get_valid_inputs() == ["matrix", "vector", "timeseries"] + + +def test_single_output(tmp_path: Path) -> None: + """Test single output setup. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_single_output.hdf5" + # Single storage, must be the uri + storage = HDF5FeatureStorage(uri=uri, single_output=True) + assert storage.single_output is True + + +def test_single_output_meta_not_found_error(tmp_path: Path) -> None: + """Test single output metadata not found error. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_single_output_no_meta.hdf5" + storage = HDF5FeatureStorage(uri=uri, single_output=True) + # Store data to create the file + storage._store_data( + kind="vector", + meta_md5="md5", + element=[{"sub": "001"}], + data=np.empty((1, 1)), + ) + # Check metadata error + with pytest.raises(IOError, match="`meta` not found in:"): + storage._read_metadata() + + +def test_single_output_file_not_found_error(tmp_path: Path) -> None: + """Test single output file not found error. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_single_output_no_file.hdf5" + storage = HDF5FeatureStorage(uri=uri, single_output=True) + # Check file error + with pytest.raises(IOError, match="HDF5 file not found at:"): + storage._read_data(md5="md5") + + +def test_multi_output_error(tmp_path: Path) -> None: + """Test error for multi output. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_multi_output.hdf5" + storage = HDF5FeatureStorage(uri=uri, single_output=False) + with pytest.raises(RuntimeError, match="`element` must be provided"): + storage._fetch_correct_uri_for_io(element=None) + + +def test_single_output_parent_path_creation(tmp_path: Path) -> None: + """Test path creation with single output creation. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + to_create_hdf5 = tmp_path / "to_create_hdf5" + # Path does not exist yet + assert not to_create_hdf5.exists() + + uri = to_create_hdf5.absolute() / "test_single_output.hdf5" + _ = HDF5FeatureStorage(uri=uri, single_output=True) + # Path exists now + assert to_create_hdf5.exists() + + +def test_store_metadata_and_list_features(tmp_path: Path) -> None: + """Test metadata store and features listing. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_metadata_store.hdf5" + # Single storage, must be the uri + storage = HDF5FeatureStorage(uri=uri, single_output=True) + # Metadata to store + meta = { + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "test"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + # List the stored features + features = storage.list_features() + # Get the first MD5 + feature_md5 = list(features.keys())[0] + # Check the MD5 + assert meta_md5 == feature_md5 + + +def test_store_metadata_ignore_duplicate(tmp_path: Path) -> None: + """Test duplicate ignore for metadata store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_duplicate_metadata_store.hdf5" + # Single storage, must be the uri + storage = HDF5FeatureStorage(uri=uri, single_output=True) + # Store metadata first time + storage.store_metadata( + meta_md5="md5", + element={"sub": "001"}, + meta={ + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "test"}, + "type": "BOLD", + }, + ) + # Store metadata second time, should be ignored + storage.store_metadata( + meta_md5="md5", + element={"sub": "001"}, + meta={ + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "test"}, + "type": "BOLD", + }, + ) + # List the stored features + features = storage.list_features() + # Should only have one element + assert len(features) == 1 + + +def test_read_df_params_error(tmp_path: Path) -> None: + """Test parameter validation errors for read_df. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_read_df_params_error.hdf5" + storage = HDF5FeatureStorage(uri=uri, single_output=True) + # Store metadata to create the file + storage.store_metadata( + meta_md5="meta_md5", + element={"sub": "001"}, + meta={ + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "test"}, + "name": "BOLD", + }, + ) + + with pytest.raises(ValueError, match="Only one of"): + storage.read_df(feature_name="name", feature_md5="md5") + + with pytest.raises(ValueError, match="At least one of"): + storage.read_df() + + with pytest.raises(ValueError, match="Feature MD5"): + storage.read_df(feature_md5="md5") + + +def test_read_df(tmp_path: Path) -> None: + """Test read_df. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_read_df.hdf5" + storage = HDF5FeatureStorage(uri=uri) + # Metadata to store + meta = { + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "test"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + # Data to store + data = np.array([[1, 10]]) + col_headers = ["f1", "f2"] + # Store table + storage.store_vector( + meta_md5=meta_md5, + element=element_to_store, + data=data, + col_names=col_headers, + ) + # Read into dataframe and check + df_md5 = storage.read_df(feature_md5=meta_md5) + df_name = storage.read_df(feature_name="BOLD_test") + assert_frame_equal(df_md5, df_name) + + # Check for errors about no / duplicate feature name + with pytest.raises(ValueError, match="not found"): + storage.read_df(feature_name="BOLD") + + # Store duplicate entry and check error + storage.store_metadata( + meta_md5=meta_md5, + element={"subject": "test-clone"}, + meta=meta_to_store, + ) + storage.store_vector( + meta_md5=meta_md5, + element={"subject": "test-clone"}, + data=data, + col_names=col_headers, + ) + + +def test_store_data_ignore_duplicate(tmp_path: Path) -> None: + """Test duplicate ignore for data store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_duplicate_data_store.hdf5" + # Single storage, must be the uri + storage = HDF5FeatureStorage(uri=uri, single_output=True) + # Store data first time + storage._store_data( + kind="vector", + meta_md5="md5", + element=[{"sub": "001"}], + data=np.empty((1, 1)), + ) + # Store data second time, should be ignored + storage._store_data( + kind="vector", + meta_md5="md5", + element=[{"sub": "001"}], + data=np.empty((1, 1)), + ) + + +def test_store_data_incorrect_kwargs(tmp_path: Path) -> None: + """Test incorrect kwargs for data store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_incorrect_kwargs_data_store.hdf5" + # Single storage, must be the uri + storage = HDF5FeatureStorage(uri=uri, single_output=True) + # Store data first time + storage._store_data( + kind="vector", + meta_md5="md5", + element=[{"sub": "001"}], + data=np.empty((1, 1)), + ) + # Store data second time, should be ignored + with pytest.raises(RuntimeError, match="The additional data for"): + storage._store_data( + kind="vector", + meta_md5="md5", + element=[{"sub": "001"}], + data=np.empty((1, 1)), + col_names="col", + ) + + +@pytest.mark.parametrize( + "force, dtype", + [ + (True, "float32"), + (False, "float64"), + ], +) +def test_f64_to_f632_conversion( + tmp_path: Path, force: bool, dtype: str +) -> None: + """Test actual data is casted from float64 to float32. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + force : bool + The parametrized conversion option. + dtype : str + The parametrized expected data type. + + """ + uri = tmp_path / "test_data_conversion.hdf5" + storage = HDF5FeatureStorage(uri=uri, force_float32=force) + # Metadata to store + meta = { + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "mark"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store matrix + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + # Store data + storage.store_matrix( + meta_md5=meta_md5, + element=element_to_store, + data=np.arange(4, dtype="float64").reshape((2, 2)), + ) + # Read into dataframe + read_df = storage.read_df(feature_md5=meta_md5) + # Check data type + assert read_df.values.dtype == np.dtype(dtype) + + +def test_store_matrix(tmp_path: Path) -> None: + """Test matrix store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_store_matrix.hdf5" + storage = HDF5FeatureStorage(uri=uri) + # Metadata to store + meta = { + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + + # Store 4 X 3 full matrix + data = np.array( + [[1, 2, 3], [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]] + ) + row_headers = ["row1", "row2", "row3", "row4"] + col_headers = ["col1", "col2", "col3"] + + # Store matrix + storage.store_matrix( + meta_md5=meta_md5, + element=element_to_store, + data=data, + row_names=row_headers, + col_names=col_headers, + ) + + # List the stored features + features = storage.list_features() + # Check the MD5 + assert "BOLD_fc" == features[meta_md5]["name"] + + # Read into dataframe + read_df = storage.read_df(feature_md5=meta_md5) + # Check shape of dataframe + assert read_df.shape == (1, 12) + # Check data of dataframe + assert_array_equal(read_df.values, data.reshape(1, -1)) + # Check column headers + assert read_df.columns.to_list() == [ + f"{row}~{col}" for row in row_headers for col in col_headers + ] + + +def test_store_matrix_without_headers(tmp_path: Path) -> None: + """Test matrix store without headers. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_store_matrix_no_headers.hdf5" + storage = HDF5FeatureStorage(uri=uri) + # Metadata to store + meta = { + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + + # Store 4 X 3 full matrix + data = np.array( + [[1, 2, 3], [11, 22, 33], [111, 222, 333], [1111, 2222, 3333]] + ) + + # Store matrix + storage.store_matrix( + meta_md5=meta_md5, element=element_to_store, data=data + ) + + # List the stored features + features = storage.list_features() + # Check the MD5 + assert "BOLD_fc" == features[meta_md5]["name"] + + # Read the dataframe + read_df = storage.read_df(feature_md5=meta_md5) + # Check shape of dataframe + assert read_df.shape == (1, 12) + # Check data of dataframe + assert_array_equal(read_df.values, data.reshape(1, -1)) + # Check column headers + assert read_df.columns.to_list() == [ + f"{row}~{col}" + for row in ["r0", "r1", "r2", "r3"] + for col in ["c0", "c1", "c2"] + ] + + +def test_store_upper_triangular_matrix(tmp_path: Path) -> None: + """Test upper triangular matrix store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_store_matrix_triu.hdf5" + storage = HDF5FeatureStorage(uri=uri) + # Metadata to store + meta = { + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + + # Store upper triangular matrix + data = np.array([[1, 2, 3], [11, 22, 33], [111, 222, 333]]) + row_headers = ["row1", "row2", "row3"] + col_headers = ["col1", "col2", "col3"] + + # Store matrix + storage.store_matrix( + meta_md5=meta_md5, + element=element_to_store, + data=data, + row_names=row_headers, + col_names=col_headers, + matrix_kind="triu", + ) + + # List the stored features + features = storage.list_features() + # Check the MD5 + assert "BOLD_fc" == features[meta_md5]["name"] + + # Read into dataframe + read_df = storage.read_df(feature_md5=meta_md5) + # Check data of dataframe + assert_array_equal(read_df.values, np.array([[1, 2, 3, 22, 33, 333]])) + # Check column headers + assert read_df.columns.to_list() == [ + "row1~col1", + "row1~col2", + "row1~col3", + "row2~col2", + "row2~col3", + "row3~col3", + ] + + +def test_store_upper_triangular_matrix_without_diagonal( + tmp_path: Path, +) -> None: + """Test upper triangular matrix without diagonal store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_store_matrix_triu_no_diagonal.hdf5" + storage = HDF5FeatureStorage(uri=uri) + # Metadata to store + meta = { + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + + # Store upper triangular matrix without diagonal + data = np.array([[1, 2, 3], [11, 22, 33], [111, 222, 333]]) + row_headers = ["row1", "row2", "row3"] + col_headers = ["col1", "col2", "col3"] + + # Store matrix + storage.store_matrix( + meta_md5=meta_md5, + element=element_to_store, + data=data, + row_names=row_headers, + col_names=col_headers, + matrix_kind="triu", + diagonal=False, + ) + + # List the stored features + features = storage.list_features() + # Check the MD5 + assert "BOLD_fc" == features[meta_md5]["name"] + + # Read into dataframe + read_df = storage.read_df(feature_md5=meta_md5) + # Check data of dataframe + assert_array_equal(read_df.values, np.array([[2, 3, 33]])) + # Check column headers + assert read_df.columns.to_list() == ["row1~col2", "row1~col3", "row2~col3"] + + +def test_store_lower_triangular_matrix(tmp_path: Path) -> None: + """Test lower triangular matrix store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_store_matrix_tril.hdf5" + storage = HDF5FeatureStorage(uri=uri) + # Metadata to store + meta = { + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + + # Store lower triangular matrix + data = np.array([[1, 2, 3], [11, 22, 33], [111, 222, 333]]) + row_headers = ["row1", "row2", "row3"] + col_headers = ["col1", "col2", "col3"] + + # Store matrix + storage.store_matrix( + meta_md5=meta_md5, + element=element_to_store, + data=data, + row_names=row_headers, + col_names=col_headers, + matrix_kind="tril", + ) + + # List the stored features + features = storage.list_features() + # Check the MD5 + assert "BOLD_fc" == features[meta_md5]["name"] + + # Read into dataframe + read_df = storage.read_df(feature_md5=meta_md5) + # Check data of dataframe + assert_array_equal(read_df.values, np.array([[1, 11, 22, 111, 222, 333]])) + # Check column headers + assert read_df.columns.to_list() == [ + "row1~col1", + "row2~col1", + "row2~col2", + "row3~col1", + "row3~col2", + "row3~col3", + ] + + +def test_store_lower_triangular_matrix_without_diagonal( + tmp_path: Path, +) -> None: + """Test lower triangular matrix without diagonal store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_store_matrix_tril_no_diagonal.hdf5" + storage = HDF5FeatureStorage(uri=uri) + # Metadata to store + meta = { + "element": {"subject": "test"}, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + + # Store lower triangular matrix without diagonal + data = np.array([[1, 2, 3], [11, 22, 33], [111, 222, 333]]) + row_headers = ["row1", "row2", "row3"] + col_headers = ["col1", "col2", "col3"] + + # Store matrix + storage.store_matrix( + meta_md5=meta_md5, + element=element_to_store, + data=data, + row_names=row_headers, + col_names=col_headers, + matrix_kind="tril", + diagonal=False, + ) + + # List the stored features + features = storage.list_features() + # Check the MD5 + assert "BOLD_fc" == features[meta_md5]["name"] + + # Read into dataframe + read_df = storage.read_df(feature_md5=meta_md5) + # Check data of dataframe + assert_array_equal(read_df.values, np.array([[11, 111, 222]])) + # Check column headers + assert read_df.columns.to_list() == ["row2~col1", "row3~col1", "row3~col2"] + + +def test_store_vector(tmp_path: Path) -> None: + """Test vector store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_store_vector.hdf5" + storage = HDF5FeatureStorage(uri=uri) + # Metadata to store + element = {"subject": "test"} + meta = { + "element": element, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + + # Data to store + data = [10, 20, 30, 40, 50] + col_names = ["f1", "f2", "f3", "f4", "f5"] + + # Store vector + storage.store_vector( + meta_md5=meta_md5, + element=element_to_store, + data=data, + col_names=col_names, + ) + + # Read into dataframe + read_df = storage.read_df(feature_md5=meta_md5) + # Check if data are equal + assert read_df.values.flatten().tolist() == data + + +def test_store_timeseries(tmp_path: Path) -> None: + """Test timeseries store. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_store_timeseries.hdf5" + storage = HDF5FeatureStorage(uri=uri) + # Metadata to store + element = {"subject": "test"} + meta = { + "element": element, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + # Process the metadata + meta_md5, meta_to_store, element_to_store = process_meta(meta) + # Store metadata + storage.store_metadata( + meta_md5=meta_md5, element=element_to_store, meta=meta_to_store + ) + + # Data to store + data = np.array([[10], [20], [30], [40], [50]]) + col_names = ["signal"] + + # Store vector + storage.store_timeseries( + meta_md5=meta_md5, + element=element_to_store, + data=data, + col_names=col_names, + ) + + # Read into dataframe + read_df = storage.read_df(feature_md5=meta_md5) + # Check if data are equal + assert_array_equal(read_df.values, data) + + +def test_multi_output_store_and_collect(tmp_path: Path): + """Test multi output storing and collection. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + uri = tmp_path / "test_multi_output_store_and_collect.hdf5" + storage = HDF5FeatureStorage(uri=uri, single_output=False) + + # Metadata to store + meta_1 = { + "element": {"subject": "test-01", "session": "ses-01"}, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + meta_2 = { + "element": {"subject": "test-02", "session": "ses-01"}, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + meta_3 = { + "element": {"subject": "test-01", "session": "ses-02"}, + "dependencies": ["numpy"], + "marker": {"name": "fc"}, + "type": "BOLD", + } + + # Data to store + data_1 = np.array([10, 20, 30, 40, 50]) + data_2 = data_1 * 10 + data_3 = data_1 * 20 + col_headers = ["f1", "f2", "f3", "f4", "f5"] + + # Process metadata for storage + hash_1, meta_to_store_1, element_to_store_1 = process_meta(meta_1) + + # Process metadata for storage + hash_2, meta_to_store_2, element_to_store_2 = process_meta(meta_2) + + # Process metadata for storage + hash_3, meta_to_store_3, element_to_store_3 = process_meta(meta_3) + + # Check hash equality as element is not considered for hash + assert hash_1 == hash_2 + assert hash_2 == hash_3 + + # Store metadata for tables + storage.store_metadata( + meta_md5=hash_1, + element=element_to_store_1, + meta=meta_to_store_1, + ) + storage.store_metadata( + meta_md5=hash_2, + element=element_to_store_2, + meta=meta_to_store_2, + ) + storage.store_metadata( + meta_md5=hash_3, + element=element_to_store_3, + meta=meta_to_store_3, + ) + + # Store tables + storage.store_vector( + meta_md5=hash_1, + element=element_to_store_1, + data=data_1, + col_names=col_headers, + ) + storage.store_vector( + meta_md5=hash_2, + element=element_to_store_2, + data=data_2, + col_names=col_headers, + ) + storage.store_vector( + meta_md5=hash_3, + element=element_to_store_3, + data=data_3, + col_names=col_headers, + ) + + # Check that base URI does not exist yet + assert not uri.exists() + + # Convert element to preifx + prefix_1 = element_to_prefix(meta_1["element"]) # type: ignore + prefix_2 = element_to_prefix(meta_2["element"]) # type: ignore + prefix_3 = element_to_prefix(meta_3["element"]) # type: ignore + + # URIs for data storage + uri_1 = uri.parent / f"{prefix_1}{uri.name}" + uri_2 = uri.parent / f"{prefix_2}{uri.name}" + uri_3 = uri.parent / f"{prefix_3}{uri.name}" + + # Check URIs for data storage exist + assert uri_1.exists() + assert uri_2.exists() + assert uri_3.exists() + + # Read stored metadata from different files using element + read_meta_1 = storage._read_metadata(element=meta_1["element"]) + read_meta_2 = storage._read_metadata(element=meta_2["element"]) + read_meta_3 = storage._read_metadata(element=meta_3["element"]) + + # Check if metadata are equal + assert read_meta_1 == read_meta_2 + assert read_meta_2 == read_meta_3 + + # Collect data + storage.collect() + # Check that base URI exists now + assert uri.exists() + + # Read unified metadata + read_unified_meta = storage.list_features() + + # Check if aggregated metadata are equal + assert read_unified_meta == {**read_meta_1, **read_meta_2, **read_meta_3} + + +def test_collect_error_single_output() -> None: + """Test error for collect in single output storage.""" + with pytest.raises( + NotImplementedError, + match="is not implemented for single output.", + ): + storage = HDF5FeatureStorage(uri="/tmp", single_output=True) + storage.collect() diff --git a/junifer/storage/tests/test_utils.py b/junifer/storage/tests/test_utils.py index 35a6f9c4b..dd7bab3bd 100644 --- a/junifer/storage/tests/test_utils.py +++ b/junifer/storage/tests/test_utils.py @@ -4,14 +4,18 @@ # Synchon Mandal # License: AGPL -from typing import Dict, List +from typing import Dict, Iterable, List, Tuple, Union +import numpy as np import pytest +from numpy.testing import assert_array_equal from junifer.storage.utils import ( element_to_prefix, get_dependency_version, + matrix_to_vector, process_meta, + store_matrix_checks, ) @@ -20,10 +24,10 @@ from junifer.storage.utils import ( [ ("click", "8.2"), ("numpy", "1.24"), - ("datalad", "0.18"), + ("datalad", "0.19"), ("pandas", "1.6"), ("nibabel", "4.1"), - ("nilearn", "1.0"), + ("nilearn", "0.10.0"), ("sqlalchemy", "1.5.0"), ("pyyaml", "7.0"), ], @@ -40,7 +44,10 @@ def test_get_dependency_version(dependency: str, max_version: str) -> None: """ version = get_dependency_version(dependency) - assert version < max_version + if len(version.split(".")) == 3: # semver + assert int(version.split(".")[1]) <= int(max_version.split(".")[1]) + else: + assert version < max_version def test_get_dependency_version_invalid() -> None: @@ -245,3 +252,171 @@ def test_element_to_prefix_invalid_type() -> None: element = 2.3 with pytest.raises(ValueError, match=r"must be a dict"): element_to_prefix(element) # type: ignore + + +@pytest.mark.parametrize( + "params, err_msg", + [ + ( + { + "matrix_kind": "half", + "diagonal": True, + "data_shape": (1, 1), + "row_names_len": 1, + "col_names_len": 1, + }, + "Invalid kind", + ), + ( + { + "matrix_kind": "full", + "diagonal": False, + "data_shape": (1, 1), + "row_names_len": 1, + "col_names_len": 1, + }, + "Diagonal cannot", + ), + ( + { + "matrix_kind": "triu", + "diagonal": False, + "data_shape": (2, 1), + "row_names_len": 2, + "col_names_len": 1, + }, + "Cannot store a non-square", + ), + ( + { + "matrix_kind": "tril", + "diagonal": False, + "data_shape": (1, 2), + "row_names_len": 1, + "col_names_len": 2, + }, + "Cannot store a non-square", + ), + ( + { + "matrix_kind": "full", + "diagonal": True, + "data_shape": (2, 2), + "row_names_len": 1, + "col_names_len": 2, + }, + "Number of row names", + ), + ( + { + "matrix_kind": "full", + "diagonal": True, + "data_shape": (2, 2), + "row_names_len": 2, + "col_names_len": 1, + }, + "Number of column names", + ), + ], +) +def test_store_matrix_checks( + params: Dict[str, Union[str, bool, Tuple[int, int], int]], err_msg: str +) -> None: + """Test matrix storing parameter checks. + + Parameters + ---------- + params : dict + The parametrized parameters for the function. + err_msg : str + The parametrized substring of expected error message. + + """ + with pytest.raises(ValueError, match=f"{err_msg}"): + store_matrix_checks(**params) # type: ignore + + +@pytest.mark.parametrize( + "params, expected_data, expected_columns", + [ + ( + { + "data": np.arange(9).reshape(3, 3), + "col_names": ["c0", "c1", "c2"], + "row_names": ("r0", "r1", "r2"), + "matrix_kind": "triu", + "diagonal": True, + }, + np.array([0, 1, 2, 4, 5, 8]), + ["r0~c0", "r0~c1", "r0~c2", "r1~c1", "r1~c2", "r2~c2"], + ), + ( + { + "data": np.arange(9).reshape(3, 3), + "col_names": ("c0", "c1", "c2"), + "row_names": ["r0", "r1", "r2"], + "matrix_kind": "triu", + "diagonal": False, + }, + np.array([1, 2, 5]), + ["r0~c1", "r0~c2", "r1~c2"], + ), + ( + { + "data": np.arange(9).reshape(3, 3), + "col_names": ("c0", "c1", "c2"), + "row_names": ["r0", "r1", "r2"], + "matrix_kind": "tril", + "diagonal": True, + }, + np.array([0, 3, 4, 6, 7, 8]), + ["r0~c0", "r1~c0", "r1~c1", "r2~c0", "r2~c1", "r2~c2"], + ), + ( + { + "data": np.arange(9).reshape(3, 3), + "col_names": ("c0", "c1", "c2"), + "row_names": ["r0", "r1", "r2"], + "matrix_kind": "tril", + "diagonal": False, + }, + np.array([3, 6, 7]), + ["r1~c0", "r2~c0", "r2~c1"], + ), + ( + { + "data": np.arange(9).reshape(3, 3), + "col_names": ("c0", "c1", "c2"), + "row_names": ["r0", "r1", "r2"], + "matrix_kind": "full", + "diagonal": False, + }, + np.arange(9), + [ + f"{r}~{c}" + for r in ["r0", "r1", "r2"] + for c in ["c0", "c1", "c2"] + ], + ), + ], +) +def test_matrix_to_vector( + params: Dict[str, Union[np.ndarray, Iterable[str], str, bool]], + expected_data: np.ndarray, + expected_columns: List[str], +) -> None: + """Test matrix to vector. + + Parameters + ---------- + params : dict + The parametrized parameters for the function. + expected_data : np.ndarray + The parametrized vector data to expect. + expected_columns : np.ndarray + The parametrized columns labels to expect. + + """ + data, columns = matrix_to_vector(**params) # type: ignore + assert_array_equal(data, expected_data) + assert columns == expected_columns diff --git a/junifer/storage/utils.py b/junifer/storage/utils.py index 05de657ba..629b64f20 100644 --- a/junifer/storage/utils.py +++ b/junifer/storage/utils.py @@ -7,7 +7,9 @@ import hashlib import json from importlib.metadata import PackageNotFoundError, version -from typing import Dict, Tuple +from typing import Dict, Iterable, List, Tuple + +import numpy as np from ..utils.logging import logger, raise_error @@ -139,3 +141,120 @@ def element_to_prefix(element: Dict) -> str: logger.debug(f"Converted prefix: {prefix}") return f"{prefix}_" + + +def store_matrix_checks( + matrix_kind: str, + diagonal: bool, + data_shape: Tuple[int, int], + row_names_len: int, + col_names_len: int, +) -> None: + """Run parameter checks for store_matrix() methods. + + Parameters + ---------- + matrix_kind : {"triu", "tril", "full"} + The kind of matrix: + + * ``triu`` : store upper triangular only + * ``tril`` : store lower triangular + * ``full`` : full matrix + + diagonal : bool + Whether to store the diagonal. If ``matrix_kind`` is "full", + setting this to False will raise an error. + data_shape : tuple of int and int + The shape of the matrix data to store. + row_names_len : int + The length of row labels. + col_names_len : int + The length of column labels. + + """ + # Matrix kind validation + if matrix_kind not in ("triu", "tril", "full"): + raise_error(msg=f"Invalid kind {matrix_kind}", klass=ValueError) + # Diagonal validation + if diagonal is False and matrix_kind not in ["triu", "tril"]: + raise_error( + msg="Diagonal cannot be False if kind is not full", + klass=ValueError, + ) + # Matrix kind and shape validation + if matrix_kind in ["triu", "tril"]: + if data_shape[0] != data_shape[1]: + raise_error( + "Cannot store a non-square matrix as a triangular matrix", + klass=ValueError, + ) + # Row label validation + if row_names_len != data_shape[0]: # type: ignore + raise_error( + msg="Number of row names does not match number of rows", + klass=ValueError, + ) + # Column label validation + if col_names_len != data_shape[1]: # type: ignore + raise_error( + msg="Number of column names does not match number of columns", + klass=ValueError, + ) + + +def matrix_to_vector( + data: np.ndarray, + col_names: Iterable[str], + row_names: Iterable[str], + matrix_kind: str, + diagonal: bool, +) -> Tuple[np.ndarray, List[str]]: + """Convert matrix to vector based on parameters. + + Parameters + ---------- + data : 2D / 3D numpy.ndarray + The matrix / tensor data to store / read. + col_names : list or tuple of str + The column labels. + row_names : list or tuple of str + The row labels. + matrix_kind : str + The kind of matrix: + + * ``triu`` : store upper triangular only + * ``tril`` : store lower triangular + * ``full`` : full matrix + + diagonal : bool + Whether to store the diagonal. + + Returns + ------- + 1D / 2D numpy.ndarray + The vector / matrix data. + list of str + The column labels. + + """ + # Prepare data indexing based on matrix kind + if matrix_kind == "triu": + k = 0 if diagonal is True else 1 + data_idx = np.triu_indices(data.shape[0], k=k) + elif matrix_kind == "tril": + k = 0 if diagonal is True else -1 + data_idx = np.tril_indices(data.shape[0], k=k) + else: # full + data_idx = ( + np.repeat(np.arange(data.shape[0]), data.shape[1]), + np.tile(np.arange(data.shape[1]), data.shape[0]), + ) + # Subset data as 1D + flat_data = data[data_idx] + # Generate flat 1D row X column names + columns = [ + f"{row_names[i]}~{col_names[j]}" # type: ignore + for i, j in zip(data_idx[0], data_idx[1]) + ] + + return flat_data, columns diff --git a/junifer/tests/test_stats.py b/junifer/tests/test_stats.py index 71b660698..ef36b1071 100644 --- a/junifer/tests/test_stats.py +++ b/junifer/tests/test_stats.py @@ -27,7 +27,7 @@ def test_get_aggfunc_by_name(name: str, params: Optional[Dict]) -> None: Parameters ---------- name : str - The paramterized name of the method name. + The parametrized name of the method name. params : dict The parametrized parameters passed to the method. diff --git a/pyproject.toml b/pyproject.toml index 635f1781e..c7c05cbef 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,13 +41,14 @@ classifiers = [ dependencies = [ "click>=8.1.3,<8.2", "numpy>=1.22,<1.24", - "datalad>=0.15.4,<0.18", + "datalad>=0.15.4,<0.19", "pandas>=1.4.0,<1.6", "nibabel>=3.2.0,<4.1", - "nilearn>=0.9.0,<1.0", + "nilearn>=0.9.0,<=0.10.0", "sqlalchemy>=1.4.27,<= 1.5.0", "pyyaml>=5.1.2,<7.0", "importlib_metadata; python_version < '3.10'", + "h5py>=3.8.0,<3.9", ] dynamic = ["version"] @@ -76,7 +77,7 @@ docs = [ ################ [tool.setuptools] -packages = ["junifer"] +packages = ["junifer", "junifer.external.h5io.h5io"] [tool.setuptools_scm] version_scheme = "guess-next-dev" @@ -86,3 +87,12 @@ write_to = "junifer/_version.py" [tool.black] line-length = 79 target-version = ["py38"] +extend-exclude = """ +( + junifer/external/h5io +) +""" + +[tool.pytest.ini_options] +minversion = "7.0" +addopts = "--ignore=junifer/external/h5io -vv" diff --git a/tox.ini b/tox.ini index 47b352ada..0b6688fa0 100644 --- a/tox.ini +++ b/tox.ini @@ -49,7 +49,7 @@ passenv = deps = pytest commands = - pytest -vv + pytest [testenv:coverage] skip_install = false @@ -57,7 +57,7 @@ deps = pytest pytest-cov commands = - pytest --cov={envsitepackagesdir}/junifer --cov-report=xml -vv --cov-report=term + pytest --cov={envsitepackagesdir}/junifer --cov-report=xml --cov-report=term [testenv:codespell] skip_install = true @@ -73,6 +73,7 @@ commands = [isort] skip = __init__.py + junifer/external/h5io profile = black line_length = 79 lines_after_imports = 2 @@ -91,6 +92,7 @@ known_third_party = [flake8] exclude = __init__.py + junifer/external/h5io max-line-length = 79 extend-ignore = ; Use of `functools.lru_cache` or `functools.cache` on methods can lead to @@ -137,6 +139,7 @@ omit = */_version.py */tests/* */junifer/configs/* + */junifer/external/h5io/* parallel = false [coverage:report] @@ -150,7 +153,7 @@ exclude_lines = precision = 2 [codespell] -skip = docs/auto_*,*.html,.git/,*.pyc,docs/_build +skip = docs/auto_*,*.html,.git/,*.pyc,docs/_build,junifer/external/h5io/ count = quiet-level = 3 ignore-words = ignore_words.txt