diff --git a/src/main/python/systemds/scuro/dataloader/tabular_loader.py b/src/main/python/systemds/scuro/dataloader/tabular_loader.py new file mode 100644 index 00000000000..6652ba2b937 --- /dev/null +++ b/src/main/python/systemds/scuro/dataloader/tabular_loader.py @@ -0,0 +1,78 @@ +# ------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +# ------------------------------------------------------------- +from dataclasses import dataclass +import numpy as np +from typing import List, Optional, Union + +from systemds.scuro.dataloader.base_loader import BaseLoader +from systemds.scuro.modality.type import ModalityType + + +@dataclass +class TabularStats: + num_instances: int + num_features: int + output_shape: tuple + output_shape_is_known: bool = True + + +class TabularLoader(BaseLoader): + def __init__( + self, + source_path: str, + indices: List[str], + feature_names: Optional[List[str]] = None, + data_type: Union[np.dtype, str] = np.float32, + chunk_size: Optional[int] = None, + normalize: bool = False, + file_format: str = "npy", + modality_type: Optional[ModalityType] = ModalityType.EMBEDDING, + ): + super().__init__(source_path, indices, data_type, chunk_size, modality_type) + self.feature_names = feature_names + self.normalize = normalize + self.file_format = file_format.lower() + if self.file_format != "npy": + raise ValueError(f"Unsupported file format: {self.file_format}") + self.stats = self.get_stats(source_path) + + def extract(self, file: str, index: Optional[Union[str, List[str]]] = None): + self.file_sanity_check(file) + data = np.load(file).astype(self._data_type, copy=False).reshape(-1) + + if self.normalize: + mean = np.mean(data) + std = np.std(data) + data = (data - mean) / (std + 1e-8) + + self.metadata.append(self.modality_type.create_metadata(data)) + self.data.append(data) + + def get_stats(self, source_path: str) -> TabularStats: + num_instances = 0 + num_features = 0 + for file_name in self.indices: + file = source_path + file_name + "." + self.file_format + self.file_sanity_check(file) + data = np.load(file) + num_features = max(num_features, int(np.prod(data.shape))) + num_instances += 1 + return TabularStats(num_instances, num_features, (num_features,)) diff --git a/src/main/python/systemds/scuro/representations/aggregate.py b/src/main/python/systemds/scuro/representations/aggregate.py index 68f42d14e98..afbd6b346ea 100644 --- a/src/main/python/systemds/scuro/representations/aggregate.py +++ b/src/main/python/systemds/scuro/representations/aggregate.py @@ -89,7 +89,7 @@ def __init__(self, aggregation_function="mean", pad_modality=True, params=None): self.aggregation_function_name = aggregation_function self.parameters = { - "aggregation_function": self._aggregation_function.keys(), + "aggregation_function": list(self._aggregation_function.keys()), } def get_current_parameters(self): diff --git a/src/main/python/systemds/scuro/representations/average.py b/src/main/python/systemds/scuro/representations/average.py index f58ba0b6802..06c435d6580 100644 --- a/src/main/python/systemds/scuro/representations/average.py +++ b/src/main/python/systemds/scuro/representations/average.py @@ -18,12 +18,12 @@ # under the License. # # ------------------------------------------------------------- -import copy from typing import List import numpy as np from systemds.scuro.modality.modality import Modality +from systemds.scuro.representations.representation import RepresentationStats from systemds.scuro.representations.utils import pad_sequences from systemds.scuro.representations.fusion import Fusion @@ -41,11 +41,38 @@ def __init__(self, params=None): self.associative = True self.commutative = True - def execute(self, modalities: List[Modality], labels=None): - data = np.asarray(copy.deepcopy(modalities[0].data), dtype=float) + def execute(self, modalities: List[Modality]): + data = np.array(modalities[0].data, dtype=np.float64) for i in range(1, len(modalities)): - data += np.asarray(modalities[i].data, dtype=float) + data += np.asarray(modalities[i].data, dtype=np.float64) data /= len(modalities) - return np.array(data) + return data + + def get_output_stats(self, input_stats_list) -> RepresentationStats: + stats_list = self._fusion_input_stats(input_stats_list) + if not stats_list: + return RepresentationStats(0, (0,)) + + num_instances = max(s.num_instances for s in stats_list) + rank = len(stats_list[0].output_shape) + if rank > 0 and all(len(s.output_shape) == rank for s in stats_list): + output_shape = tuple( + max(s.output_shape[d] for s in stats_list) for d in range(rank) + ) + else: + output_shape = max(stats_list, key=self._stats_num_elements).output_shape + output_shape_is_known = all(s.output_shape_is_known for s in stats_list) + return RepresentationStats(num_instances, output_shape, output_shape_is_known) + + def estimate_peak_memory_bytes(self, input_stats) -> dict: + stats_list = self._as_stats_list(input_stats) + input_bytes = sum(self._stats_bytes(s) for s in stats_list) + output_bytes = self._stats_bytes(self.get_output_stats(input_stats)) + + raw_bytes = self._raw_input_bytes(input_stats) + cpu_peak = ( + int((raw_bytes + input_bytes + 2 * output_bytes) * 1.15) + 8 * 1024 * 1024 + ) + return {"cpu_peak_bytes": cpu_peak, "gpu_peak_bytes": 0} diff --git a/src/main/python/systemds/scuro/representations/bert.py b/src/main/python/systemds/scuro/representations/bert.py index 245466afb43..9e60f843416 100644 --- a/src/main/python/systemds/scuro/representations/bert.py +++ b/src/main/python/systemds/scuro/representations/bert.py @@ -102,6 +102,7 @@ def get_output_stats(self, input_stats) -> RepresentationStats: input_stats.num_instances, (self.max_seq_length, 768), aggregate_dim=(0,), + dtype=self.data_type, ) else: self.stats = RepresentationStats( @@ -111,6 +112,7 @@ def get_output_stats(self, input_stats) -> RepresentationStats: 0, 1, ), + dtype=self.data_type, ) if self.params and "_pushdown_aggregation" in self.params: output_shape = (768,) diff --git a/src/main/python/systemds/scuro/representations/bow.py b/src/main/python/systemds/scuro/representations/bow.py index 9e55add5de0..9a7766106c1 100644 --- a/src/main/python/systemds/scuro/representations/bow.py +++ b/src/main/python/systemds/scuro/representations/bow.py @@ -18,6 +18,8 @@ # under the License. # # ------------------------------------------------------------- +import os + import numpy as np from sklearn.feature_extraction.text import CountVectorizer @@ -30,6 +32,8 @@ from systemds.scuro.drsearch.operator_registry import register_representation from systemds.scuro.dataloader.text_loader import TextStats +_MAX_VOCAB_FEATURES = int(os.environ.get("SCURO_BOW_MAX_FEATURES", "100000")) + @register_representation(ModalityType.TEXT) class BoW(UnimodalRepresentation): @@ -43,14 +47,17 @@ def __init__(self, ngram_range=2, min_df=2, output_file=None, params=None): def get_output_stats(self, input_stats: TextStats) -> RepresentationStats: vocab_estimate = min( - 100000, + _MAX_VOCAB_FEATURES, max( 1000, input_stats.num_instances * input_stats.max_length * self.ngram_range, ), ) return RepresentationStats( - input_stats.num_instances, (vocab_estimate,), output_shape_is_known=False + input_stats.num_instances, + (vocab_estimate,), + output_shape_is_known=False, + dtype=self.data_type, ) def estimate_output_memory_bytes(self, input_stats: TextStats) -> int: @@ -73,7 +80,10 @@ def estimate_peak_memory_bytes(self, input_stats: TextStats) -> dict: def transform(self, modality, aggregation=None): transformed_modality = TransformedModality(modality, self) vectorizer = CountVectorizer( - ngram_range=(1, self.ngram_range), min_df=self.min_df + ngram_range=(1, self.ngram_range), + min_df=self.min_df, + max_features=_MAX_VOCAB_FEATURES, + dtype=np.float32, ) X = ( diff --git a/src/main/python/systemds/scuro/representations/clip.py b/src/main/python/systemds/scuro/representations/clip.py index 2c880686b11..c4e28404466 100644 --- a/src/main/python/systemds/scuro/representations/clip.py +++ b/src/main/python/systemds/scuro/representations/clip.py @@ -109,13 +109,17 @@ def get_output_stats(self, input_stats) -> RepresentationStats: input_stats.max_length, 512, ), + dtype=self.data_type, ) elif not isinstance(input_stats, RepresentationStats): - return RepresentationStats(input_stats.num_instances, (512,)) + return RepresentationStats( + input_stats.num_instances, (512,), dtype=self.data_type + ) else: return RepresentationStats( input_stats.num_instances, (input_stats.output_shape[0], 512), + dtype=self.data_type, ) def estimate_peak_memory_bytes(self, input_stats) -> dict: @@ -394,7 +398,10 @@ def _get_parameters(self): def get_output_stats(self, input_stats) -> RepresentationStats: if not isinstance(input_stats, RepresentationStats): self.stats = RepresentationStats( - input_stats.num_instances, (512,), aggregate_dim=(0,) + input_stats.num_instances, + (512,), + aggregate_dim=(0,), + dtype=self.data_type, ) else: self.stats = RepresentationStats( @@ -404,6 +411,7 @@ def get_output_stats(self, input_stats) -> RepresentationStats: 0, 1, ), + dtype=self.data_type, ) if self.params and "_pushdown_aggregation" in self.params: output_shape = (512,) diff --git a/src/main/python/systemds/scuro/representations/concatenation.py b/src/main/python/systemds/scuro/representations/concatenation.py index 3bdfdb28b1f..cee30382963 100644 --- a/src/main/python/systemds/scuro/representations/concatenation.py +++ b/src/main/python/systemds/scuro/representations/concatenation.py @@ -34,11 +34,70 @@ @register_fusion_operator() class Concatenation(Fusion): - def __init__(self, params=None): - """ - Combines modalities using concatenation - """ + def __init__(self, params=None, preserve_leading_axis=False): super().__init__("Concatenation") + if params is not None: + preserve_leading_axis = params.get( + "preserve_leading_axis", preserve_leading_axis + ) + + self.preserve_leading_axis = bool(preserve_leading_axis) + self.preserves_leading_axis = self.preserve_leading_axis + + def get_current_parameters(self): + current_params = super().get_current_parameters() + current_params["preserve_leading_axis"] = self.preserve_leading_axis + return current_params + + @staticmethod + def _as_dense(modality): + dtype = modality.metadata[0]["data_layout"]["type"] + data = modality.data + arr = ( + np.asarray(data, dtype=dtype) if not isinstance(data, np.ndarray) else data + ) + if arr.dtype == object: + instances = [np.asarray(instance, dtype=dtype) for instance in data] + rest = tuple( + max(i.shape[d] for i in instances) for d in range(instances[0].ndim) + ) + arr = np.zeros((len(instances), *rest), dtype=dtype) + for i, instance in enumerate(instances): + arr[(i, *(slice(0, s) for s in instance.shape))] = instance + return arr + + @staticmethod + def _to_window_feature_matrix(arr): + if arr.ndim == 1: + return arr[:, None, None] + if arr.ndim == 2: + return arr[:, :, None] + return arr.reshape(arr.shape[0], arr.shape[1], -1) + + @staticmethod + def _flatten_feature_shape(shape): + if len(shape) == 0: + return (1, 1) + if len(shape) == 1: + return (shape[0], 1) + return (shape[0], int(np.prod(shape[1:]))) + + def _concat_on_leading_axis(self, modalities: List[Modality]): + arrays = [ + self._to_window_feature_matrix(self._as_dense(modality)) + for modality in modalities + ] + + num_windows = max(arr.shape[1] for arr in arrays) + aligned = [] + for arr in arrays: + if arr.shape[1] < num_windows: + pad_width = [(0, 0)] * arr.ndim + pad_width[1] = (0, num_windows - arr.shape[1]) + arr = np.pad(arr, pad_width=pad_width, mode="constant") + aligned.append(arr) + + return np.concatenate(aligned, axis=-1) def execute(self, modalities: List[Modality]): if len(modalities) == 1: @@ -47,6 +106,9 @@ def execute(self, modalities: List[Modality]): dtype=modalities[0].metadata[0]["data_layout"]["type"], ) + if self.preserve_leading_axis: + return self._concat_on_leading_axis(modalities) + max_emb_size = self.get_max_embedding_size(modalities) size = len(modalities[0].data) @@ -71,35 +133,32 @@ def execute(self, modalities: List[Modality]): return np.array(data) def get_output_stats(self, input_stats_list) -> RepresentationStats: - if isinstance(input_stats_list, RepresentationStats): - return input_stats_list - - stats_list = list(input_stats_list) + stats_list = self._fusion_input_stats(input_stats_list) if not stats_list: return RepresentationStats(0, (0,)) - num_instances = stats_list[0].num_instances - total_dim = sum(s.output_shape[-1] for s in stats_list) - output_shape = (total_dim,) - - return RepresentationStats(num_instances, output_shape) + num_instances = max(s.num_instances for s in stats_list) + shapes = [tuple(s.output_shape) for s in stats_list] + if self.preserve_leading_axis: + shapes = [self._flatten_feature_shape(shape) for shape in shapes] + rank = len(shapes[0]) - def estimate_peak_memory_bytes(self, input_stats_list) -> dict: - elem_size = np.dtype(np.float32).itemsize - - def stats_bytes(s: RepresentationStats) -> int: - numel = int(np.prod(s.output_shape)) if len(s.output_shape) > 0 else 1 - return int(s.num_instances * numel * elem_size) + if rank >= 1 and all(len(shape) == rank for shape in shapes): + leading = tuple(max(shape[d] for shape in shapes) for d in range(rank - 1)) + output_shape = (*leading, sum(shape[-1] for shape in shapes)) + else: + output_shape = max(stats_list, key=self._stats_num_elements).output_shape - current_output = 0 - peak = 0 - for s in input_stats_list: - chunk = stats_bytes(s) - new_output = current_output + chunk + output_shape_is_known = all(s.output_shape_is_known for s in stats_list) + return RepresentationStats(num_instances, output_shape, output_shape_is_known) - step_peak = current_output + chunk + new_output + chunk - peak = max(peak, step_peak) - current_output = new_output + def estimate_peak_memory_bytes(self, input_stats) -> dict: + stats_list = self._as_stats_list(input_stats) + input_bytes = sum(self._stats_bytes(s) for s in stats_list) + output_bytes = self._stats_bytes(self.get_output_stats(input_stats)) - cpu_peak = int(peak * 1.15 + 16 * 1024 * 1024) + raw_bytes = self._raw_input_bytes(input_stats) + cpu_peak = ( + int((raw_bytes + 2 * input_bytes + output_bytes) * 1.1) + 8 * 1024 * 1024 + ) return {"cpu_peak_bytes": cpu_peak, "gpu_peak_bytes": 0} diff --git a/src/main/python/systemds/scuro/representations/fusion.py b/src/main/python/systemds/scuro/representations/fusion.py index 1426797f00b..e56eb465bc6 100644 --- a/src/main/python/systemds/scuro/representations/fusion.py +++ b/src/main/python/systemds/scuro/representations/fusion.py @@ -30,7 +30,14 @@ from systemds.scuro.modality.transformed import TransformedModality from systemds.scuro.modality.modality import Modality -from systemds.scuro.representations.representation import Representation +from systemds.scuro.representations.representation import ( + CONTAINER_ARRAY, + Representation, + RepresentationStats, + derive_stats, + stats_bytes, + stats_num_elements, +) from systemds.scuro.utils.schema_helpers import get_shape @@ -46,7 +53,42 @@ def __init__(self, name, parameters=None): self.needs_alignment = False self.needs_training = False self.needs_instance_alignment = False + self.preserves_leading_axis = False self.output_modality_type = ModalityType.EMBEDDING + self.data_type = np.float32 + + @staticmethod + def _as_stats_list(input_stats) -> List[RepresentationStats]: + if isinstance(input_stats, RepresentationStats): + return [input_stats] + return list(input_stats) + + def _pre_aggregated_stats(self, stats: RepresentationStats) -> RepresentationStats: + shape = tuple(stats.output_shape) + if len(shape) > 1 and not self.preserves_leading_axis: + shape = shape[1:] + + return derive_stats( + stats, + output_shape=shape, + dtype=self.data_type, + container=CONTAINER_ARRAY, + ) + + def _fusion_input_stats(self, input_stats) -> List[RepresentationStats]: + return [self._pre_aggregated_stats(s) for s in self._as_stats_list(input_stats)] + + def _raw_input_bytes(self, input_stats) -> int: + return sum( + stats_bytes(s, quantile=0.95) for s in self._as_stats_list(input_stats) + ) + + @staticmethod + def _stats_num_elements(stats: RepresentationStats) -> int: + return stats_num_elements(stats) + + def _stats_bytes(self, stats: RepresentationStats) -> int: + return stats_bytes(stats) def transform(self, modalities: List[Modality]): """ @@ -58,12 +100,14 @@ def transform(self, modalities: List[Modality]): mods = [] for modality in modalities: agg_modality = None - if get_shape(modality.metadata) > 1: + if not self.preserves_leading_axis and get_shape(modality.metadata) > 1: agg_operator = AggregatedRepresentation() agg_modality = agg_operator.transform(modality) mods.append(agg_modality if agg_modality else modality) if self.needs_alignment: + for modality in mods: + self._normalize_for_fusion(modality) max_len = self.get_max_embedding_size(mods) for modality in mods: modality.pad(max_len=max_len) @@ -124,6 +168,20 @@ def apply_representation(self, modalities: List[Modality]): else: return self.execute(modalities) + @staticmethod + def _normalize_for_fusion(modality: Modality): + if not modality.has_data(): + return + + arr = np.asarray(modality.data) + if arr.dtype == object or arr.ndim != 1: + return + + if modality.has_metadata() and len(modality.metadata) == 1: + modality.data = arr.reshape(1, -1).copy() + elif not modality.has_metadata() or len(modality.metadata) <= 1: + modality.data = arr.reshape(1, -1).copy() + def get_max_embedding_size(self, modalities: List[Modality]): """ Computes the maximum embedding size from a given list of modalities @@ -137,9 +195,19 @@ def get_max_embedding_size(self, modalities: List[Modality]): if isinstance(data, memoryview): data = np.array(data) arr = np.asarray(data) - if arr.ndim < 2: + if arr.dtype == object: + continue + if arr.ndim == 1: + if m.has_metadata() and len(m.metadata) == 1: + emb_size = arr.shape[0] + elif not m.has_metadata() or len(m.metadata) <= 1: + emb_size = arr.shape[0] + else: + continue + elif arr.ndim >= 2: + emb_size = arr.shape[1] + else: continue - emb_size = arr.shape[1] if emb_size > max_size: max_size = emb_size return max_size diff --git a/src/main/python/systemds/scuro/representations/glove.py b/src/main/python/systemds/scuro/representations/glove.py index b45213ae19a..acf7aff522b 100644 --- a/src/main/python/systemds/scuro/representations/glove.py +++ b/src/main/python/systemds/scuro/representations/glove.py @@ -59,7 +59,9 @@ def __init__(self, output_file=None, params=None): self.embedding_dim = 100 def get_output_stats(self, input_stats: TextStats) -> RepresentationStats: - return RepresentationStats(input_stats.num_instances, (self.embedding_dim,)) + return RepresentationStats( + input_stats.num_instances, (self.embedding_dim,), dtype=self.data_type + ) def estimate_output_memory_bytes(self, input_stats: TextStats) -> int: output_bytes = 1 diff --git a/src/main/python/systemds/scuro/representations/hadamard.py b/src/main/python/systemds/scuro/representations/hadamard.py index fc053f9c6dc..0697a9f67ea 100644 --- a/src/main/python/systemds/scuro/representations/hadamard.py +++ b/src/main/python/systemds/scuro/representations/hadamard.py @@ -48,29 +48,28 @@ def execute(self, modalities: List[Modality], train_indices=None): return fused_data def get_output_stats(self, input_stats_list) -> RepresentationStats: - if isinstance(input_stats_list, RepresentationStats): - return input_stats_list - - stats_list = list(input_stats_list) + stats_list = self._fusion_input_stats(input_stats_list) if not stats_list: return RepresentationStats(0, (0,)) - max_dim = max([stats.output_shape[-1] for stats in stats_list]) - return RepresentationStats(stats_list[0].num_instances, (max_dim,)) - - def estimate_peak_memory_bytes(self, input_stats_list) -> dict: - elem_size = np.dtype(np.float64).itemsize + num_instances = max(s.num_instances for s in stats_list) + rank = len(stats_list[0].output_shape) + if rank > 0 and all(len(s.output_shape) == rank for s in stats_list): + output_shape = tuple( + max(s.output_shape[d] for s in stats_list) for d in range(rank) + ) + else: + output_shape = max(stats_list, key=self._stats_num_elements).output_shape + output_shape_is_known = all(s.output_shape_is_known for s in stats_list) + return RepresentationStats(num_instances, output_shape, output_shape_is_known) - def stats_payload_bytes(s: RepresentationStats) -> int: - numel = int(np.prod(s.output_shape)) if len(s.output_shape) > 0 else 1 - return int(s.num_instances * numel * elem_size) + def estimate_peak_memory_bytes(self, input_stats) -> dict: + stats_list = self._as_stats_list(input_stats) + input_bytes = sum(self._stats_bytes(s) for s in stats_list) + output_bytes = self._stats_bytes(self.get_output_stats(input_stats)) - stacked_input_bytes = sum(stats_payload_bytes(s) for s in input_stats_list) - out_stats = self.get_output_stats(input_stats_list) - output_bytes = stats_payload_bytes(out_stats) - reduction_workspace_bytes = output_bytes - cpu_peak = int( - (stacked_input_bytes + output_bytes + reduction_workspace_bytes) * 1.15 - + 8 * 1024 * 1024 + raw_bytes = self._raw_input_bytes(input_stats) + cpu_peak = ( + int((raw_bytes + 2 * input_bytes + output_bytes) * 1.15) + 8 * 1024 * 1024 ) return {"cpu_peak_bytes": cpu_peak, "gpu_peak_bytes": 0} diff --git a/src/main/python/systemds/scuro/representations/mfcc.py b/src/main/python/systemds/scuro/representations/mfcc.py index 483ae3eef81..804dc04f1f2 100644 --- a/src/main/python/systemds/scuro/representations/mfcc.py +++ b/src/main/python/systemds/scuro/representations/mfcc.py @@ -58,7 +58,7 @@ def __init__( "n_fft": [1024, 2048, 4096], } - super().__init__("MFCC", ModalityType.TIMESERIES, parameters, False) + super().__init__("MFCC", ModalityType.TIMESERIES, parameters, True) if params is not None: n_mfcc = params.get("n_mfcc", n_mfcc) diff --git a/src/main/python/systemds/scuro/representations/mlp_averaging.py b/src/main/python/systemds/scuro/representations/mlp_averaging.py index fb71424d738..1150a149ccd 100644 --- a/src/main/python/systemds/scuro/representations/mlp_averaging.py +++ b/src/main/python/systemds/scuro/representations/mlp_averaging.py @@ -24,13 +24,11 @@ from torch.utils.data import DataLoader, TensorDataset import numpy as np -import warnings +import logging from systemds.scuro.modality.type import ModalityType from systemds.scuro.representations.representation import RepresentationStats from systemds.scuro.utils.static_variables import ( - compute_batch_size, get_device, - get_device_for_model, ) from systemds.scuro.utils.utils import set_random_seeds from systemds.scuro.drsearch.operator_registry import ( @@ -40,6 +38,8 @@ DimensionalityReduction, ) +logger = logging.getLogger(__name__) + @register_dimensionality_reduction_operator(ModalityType.EMBEDDING) class MLPAveraging(DimensionalityReduction): @@ -144,8 +144,14 @@ def execute(self, data): input_dim = data.shape[1] if input_dim <= self.output_dim: - warnings.warn( - f"Input dimension {input_dim} is smaller than output dimension {self.output_dim}. Returning original data." + # Expected outcome, not a defect: the search offers MLPAveraging + # every output_dim in its parameter grid, so a narrow input hits + # this on most of them. A warning per call buried the run log, so + # it is reported at debug level instead. + logger.debug( + "Input dimension %d is smaller than output dimension %d. Returning original data.", + input_dim, + self.output_dim, ) # TODO: this should be pruned as possible representation, could add output_dim as parameter to reps if possible return data diff --git a/src/main/python/systemds/scuro/representations/physiological_representations.py b/src/main/python/systemds/scuro/representations/physiological_representations.py new file mode 100644 index 00000000000..669732855fa --- /dev/null +++ b/src/main/python/systemds/scuro/representations/physiological_representations.py @@ -0,0 +1,641 @@ +# ------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +# ------------------------------------------------------------- +import numpy as np +from scipy.signal import find_peaks +from scipy.spatial.distance import pdist + +from systemds.scuro.modality.type import ModalityType +from systemds.scuro.modality.transformed import TransformedModality +from systemds.scuro.representations.representation import RepresentationStats +from systemds.scuro.representations.unimodal import UnimodalRepresentation +from systemds.scuro.drsearch.operator_registry import ( + register_representation, + register_context_representation_operator, +) + + +def _as_float1d(signal): + return np.asarray(signal, dtype=np.float64).reshape(-1) + + +def _successive_diffs(intervals): + intervals = _as_float1d(intervals) + if intervals.size < 2: + return np.array([], dtype=np.float64) + return np.diff(intervals) + + +def _segment_duration(intervals): + intervals = _as_float1d(intervals) + if intervals.size == 0: + return 0.0 + return float(np.sum(intervals)) + + +def _interpolate_tachogram(nn_intervals, fs=4.0): + nn = _as_float1d(nn_intervals) + if nn.size < 2: + return nn + times = np.cumsum(nn) + times = times - times[0] + duration = times[-1] + if duration <= 0.0: + return nn + t_uniform = np.arange(0.0, duration, 1.0 / fs) + if t_uniform.size < 2: + return nn + return np.interp(t_uniform, times, nn) + + +def _bandpower(signal, fs, f1, f2): + x = _as_float1d(signal) + if x.size < 4: + return 0.0 + x = x - np.mean(x) + freqs = np.fft.rfftfreq(x.size, d=1.0 / fs) + psd = np.abs(np.fft.rfft(x)) ** 2 + mask = (freqs >= f1) & (freqs < f2) + return float(np.sum(psd[mask])) + + +def _detect_ecg_r_peaks(signal, fs, min_rr_s=0.3): + x = _as_float1d(signal) + min_distance = max(1, int(min_rr_s * fs)) + + if x.size <= min_distance: + return np.array([], dtype=int) + # 99th percentile instead of max to stay robust to outlier spikes. + height = np.quantile(x, 0.99) / 2.0 + peaks, _ = find_peaks(x, distance=min_distance, height=height) + return peaks + + +def _ecg_nn_intervals(signal, fs): + peaks = _detect_ecg_r_peaks(signal, fs) + if peaks.size < 2: + return np.array([], dtype=np.float64) + return np.diff(peaks) / float(fs) + + +def _count_matches_vectorized(x, template_len, r): + n = x.size + n_templates = n - template_len + if n_templates < 2: + return 0 + templates = np.lib.stride_tricks.sliding_window_view(x, template_len)[:n_templates] + dists = pdist(templates, metric="chebyshev") + return int(np.sum(dists <= r)) + + +def _sample_entropy(intervals, m=2, r_factor=0.2): + x = _as_float1d(intervals) + n = x.size + if n <= m + 1: + return 0.0 + r = r_factor * np.std(x) + if r <= 0.0: + return 0.0 + + a = _count_matches_vectorized(x, m + 1, r) + b = _count_matches_vectorized(x, m, r) + if b == 0 or a == 0: + return 0.0 + return float(-np.log(a / b)) + + +def _fluctuation_for_scale(y, scale, n_segments): + segments = y[: n_segments * scale].reshape(n_segments, scale) + t = np.arange(scale, dtype=np.float64) + t_mean = t.mean() + t_centered = t - t_mean + denom = np.sum(t_centered**2) + if denom == 0.0: + return None + seg_mean = segments.mean(axis=1, keepdims=True) + slope = (segments * t_centered[None, :]).sum(axis=1, keepdims=True) / denom + intercept = seg_mean - slope * t_mean + trend = intercept + slope * t[None, :] + rms = np.sqrt(np.mean((segments - trend) ** 2, axis=1)) + return float(rms.mean()) + + +def _dfa_alpha(signal, min_box=4, max_box=None): + x = _as_float1d(signal) + n = x.size + if n < min_box * 2: + return 0.0 + y = np.cumsum(x - np.mean(x)) + if max_box is None: + max_box = max(min_box + 1, n // 4) + if max_box <= min_box: + return 0.0 + + scales = np.unique( + np.logspace(np.log10(min_box), np.log10(max_box), num=10, dtype=int) + ) + fluctuations = [] + for scale in scales: + if scale < min_box: + continue + n_segments = n // scale + if n_segments < 1: + continue + fluctuation = _fluctuation_for_scale(y, int(scale), n_segments) + if fluctuation is not None: + fluctuations.append((scale, fluctuation)) + + fluctuations = [(s, f) for s, f in fluctuations if f > 0.0] + if len(fluctuations) < 2: + return 0.0 + scales, f_vals = zip(*fluctuations) + alpha = np.polyfit(np.log(scales), np.log(f_vals), 1)[0] + return float(alpha) + + +def _detect_scr_peaks(signal, fs, min_distance_s=1.0, prominence_factor=0.05): + x = _as_float1d(signal) + min_distance = max(1, int(min_distance_s * fs)) + + if x.size <= min_distance: + return ( + np.array([], dtype=int), + np.array([], dtype=float), + np.array([], dtype=float), + ) + + prominence = max(prominence_factor * (np.max(x) - np.min(x)), 1e-12) + peaks, props = find_peaks(x, distance=min_distance, prominence=prominence) + + if peaks.size == 0: + return peaks, np.array([], dtype=np.float64), np.array([], dtype=np.float64) + + left = np.maximum(0, peaks - min_distance) + right = np.minimum(x.size - 1, peaks + min_distance) + window = min_distance + 1 + offsets = np.arange(window) + + k_max_left = peaks - left + idx_left = np.clip(peaks[:, None] - offsets[None, :], 0, x.size - 1) + vals_left = x[idx_left] + valid_left = offsets[None, :] <= k_max_left[:, None] + + baseline = np.min(np.where(valid_left, vals_left, np.inf), axis=1) + amplitudes = x[peaks] - baseline + half = baseline + 0.5 * amplitudes + + below_left = (vals_left <= half[:, None]) & valid_left + has_below_left = below_left.any(axis=1) + first_below_left = np.argmax(below_left, axis=1) + li = np.maximum( + left, peaks - np.where(has_below_left, first_below_left, k_max_left) + ) + + m_max = right - peaks + idx_right = np.clip(peaks[:, None] + offsets[None, :], 0, x.size - 1) + vals_right = x[idx_right] + valid_right = offsets[None, :] <= m_max[:, None] + + below_right = (vals_right <= half[:, None]) & valid_right + has_below_right = below_right.any(axis=1) + first_below_right = np.argmax(below_right, axis=1) + ri = np.minimum(right, peaks + np.where(has_below_right, first_below_right, m_max)) + + durations = (ri - li) / fs + + return ( + peaks, + amplitudes.astype(np.float64), + durations.astype(np.float64), + ) + + +class PhysiologicalRepresentation(UnimodalRepresentation): + def __init__(self, name, parameters=None, params=None): + if params is None: + params = {} + super().__init__(name, ModalityType.EMBEDDING, parameters, False) + + def compute_feature(self, signal): + raise NotImplementedError("Subclasses should implement this method.") + + def transform(self, modality, aggregation=None): + transformed_modality = TransformedModality( + modality, self, self.output_modality_type + ) + result = [] + for signal in modality.data: + result.append(self.compute_feature(signal)) + + maxlen = max(r.size for r in result) + padded_result = [ + np.pad(r, (0, maxlen - r.size), mode="constant", constant_values=0.0) + for r in result + ] + dtype = modality.metadata[0]["data_layout"]["type"] + transformed_modality.data = np.vstack(np.asarray(padded_result)).astype(dtype) + return transformed_modality + + def get_output_stats(self, input_stats): + return RepresentationStats( + input_stats.num_instances, (1,), input_stats.output_shape_is_known + ) + + @staticmethod + def _num_elements(shape) -> int: + n = 1 + for d in shape: + n *= int(d) + return n + + def estimate_output_memory_bytes(self, input_stats): + out_stats = self.get_output_stats(input_stats) + return ( + int(out_stats.num_instances) + * self._num_elements(out_stats.output_shape) + * np.dtype(np.float32).itemsize + ) + + def estimate_peak_memory_bytes(self, input_stats): + n_per_instance = self._num_elements(input_stats.output_shape) + input_bytes = ( + int(input_stats.num_instances) + * n_per_instance + * np.dtype(np.float32).itemsize + ) + transient_bytes = n_per_instance * np.dtype(np.float64).itemsize + output_bytes = self.estimate_output_memory_bytes(input_stats) + cpu_peak = ( + int((input_bytes + 2 * transient_bytes + output_bytes) * 1.15) + + 4 * 1024 * 1024 + ) + return {"cpu_peak_bytes": cpu_peak, "gpu_peak_bytes": 0} + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class SDNN(PhysiologicalRepresentation): + def __init__(self, fs=500.0, params=None): + super().__init__("SDNN") + if params is not None: + fs = params.get("fs", fs) + self.fs = fs + + def compute_feature(self, signal): + nn = _ecg_nn_intervals(signal, self.fs) + if nn.size < 2: + return np.array(0.0) + return np.array(np.std(nn, ddof=1)) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class RMSSD(PhysiologicalRepresentation): + def __init__(self, fs=500.0, params=None): + super().__init__("RMSSD") + if params is not None: + fs = params.get("fs", fs) + self.fs = fs + + def compute_feature(self, signal): + nn = _ecg_nn_intervals(signal, self.fs) + diffs = _successive_diffs(nn) + if diffs.size == 0: + return np.array(0.0) + return np.array(np.sqrt(np.mean(diffs**2))) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class pNN(PhysiologicalRepresentation): + def __init__(self, threshold_ms=50, fs=500.0, params=None): + super().__init__("pNN", parameters={"threshold_ms": [20, 50]}) + if params is not None: + threshold_ms = params.get("threshold_ms", threshold_ms) + fs = params.get("fs", fs) + self.threshold_ms = threshold_ms + self.fs = fs + + def compute_feature(self, signal): + nn = _ecg_nn_intervals(signal, self.fs) + diffs = _successive_diffs(nn) + if diffs.size == 0: + return np.array(0.0) + threshold = self.threshold_ms / 1000.0 + return np.array(100.0 * np.mean(np.abs(diffs) > threshold)) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class RRPerMinute(PhysiologicalRepresentation): + def __init__(self, fs=500.0, params=None): + super().__init__("RRPerMinute") + if params is not None: + fs = params.get("fs", fs) + self.fs = fs + + def compute_feature(self, signal): + nn = _ecg_nn_intervals(signal, self.fs) + duration = _segment_duration(nn) + if duration <= 0.0 or nn.size == 0: + return np.array(0.0) + return np.array(60.0 * nn.size / duration) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class HRVBandPower(PhysiologicalRepresentation): + def __init__(self, fs=4.0, f1=0.04, f2=0.15, signal_fs=500.0, params=None): + super().__init__( + "HRVBandPower", + parameters={ + "fs": [2.0, 4.0], + "f1": [0.003, 0.04, 0.15], + "f2": [0.04, 0.15, 0.4], + }, + ) + if params is not None: + fs = params.get("fs", fs) + f1 = params.get("f1", f1) + f2 = params.get("f2", f2) + signal_fs = params.get("signal_fs", signal_fs) + self.fs = fs + self.f1 = f1 + self.f2 = f2 + self.signal_fs = signal_fs + + def compute_feature(self, signal): + nn = _ecg_nn_intervals(signal, self.signal_fs) + tach = _interpolate_tachogram(nn, fs=self.fs) + return np.array(_bandpower(tach, self.fs, self.f1, self.f2)) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class HRVVLF(HRVBandPower): + def __init__(self, fs=4.0, signal_fs=500.0, params=None): + super().__init__(fs=fs, f1=0.003, f2=0.04, signal_fs=signal_fs, params=params) + self.name = "HRVVLF" + self._parameters = {"fs": [2.0, 4.0]} + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class HRVLF(HRVBandPower): + def __init__(self, fs=4.0, signal_fs=500.0, params=None): + super().__init__(fs=fs, f1=0.04, f2=0.15, signal_fs=signal_fs, params=params) + self.name = "HRVLF" + self._parameters = {"fs": [2.0, 4.0]} + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class HRVHF(HRVBandPower): + def __init__(self, fs=4.0, signal_fs=500.0, params=None): + super().__init__(fs=fs, f1=0.15, f2=0.40, signal_fs=signal_fs, params=params) + self.name = "HRVHF" + self._parameters = {"fs": [2.0, 4.0]} + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class HRVLFHF(PhysiologicalRepresentation): + def __init__(self, fs=4.0, signal_fs=500.0, params=None): + super().__init__("HRVLFHF", parameters={"fs": [2.0, 4.0]}) + if params is not None: + fs = params.get("fs", fs) + signal_fs = params.get("signal_fs", signal_fs) + self.fs = fs + self.signal_fs = signal_fs + + def compute_feature(self, signal): + lf = HRVLF(fs=self.fs, signal_fs=self.signal_fs).compute_feature(signal)[()] + hf = HRVHF(fs=self.fs, signal_fs=self.signal_fs).compute_feature(signal)[()] + if hf <= 0.0: + return np.array(0.0) + return np.array(lf / hf) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class PoincareSD1(PhysiologicalRepresentation): + def __init__(self, fs=500.0, params=None): + super().__init__("PoincareSD1") + if params is not None: + fs = params.get("fs", fs) + self.fs = fs + + def compute_feature(self, signal): + nn = _ecg_nn_intervals(signal, self.fs) + diffs = _successive_diffs(nn) + if diffs.size < 2: + return np.array(0.0) + return np.array(np.std(diffs, ddof=1) / np.sqrt(2.0)) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class PoincareSD2(PhysiologicalRepresentation): + def __init__(self, fs=500.0, params=None): + super().__init__("PoincareSD2") + if params is not None: + fs = params.get("fs", fs) + self.fs = fs + + def compute_feature(self, signal): + nn = _ecg_nn_intervals(signal, self.fs) + if nn.size < 3: + return np.array(0.0) + summed = nn[:-1] + nn[1:] + return np.array(np.std(summed, ddof=1) / np.sqrt(2.0)) + + +# @register_representation([ModalityType.PHYSIOLOGICAL]) +# @register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +# class SampleEntropy(PhysiologicalRepresentation): +# def __init__(self, m=2, r_factor=0.2, params=None): +# super().__init__( +# "SampleEntropy", parameters={"m": [2, 3], "r_factor": [0.15, 0.2, 0.25]} +# ) +# if params is not None: +# m = params.get("m", m) +# r_factor = params.get("r_factor", r_factor) +# self.m = m +# self.r_factor = r_factor + +# def compute_feature(self, nn): +# return np.array(_sample_entropy(nn, m=self.m, r_factor=self.r_factor)) + + +# @register_representation([ModalityType.PHYSIOLOGICAL]) +# @register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +# class DFAAlpha(PhysiologicalRepresentation): +# def __init__(self, params=None): +# super().__init__("DFAAlpha") + +# def compute_feature(self, nn): +# return np.array(_dfa_alpha(nn)) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class SCLSlope(PhysiologicalRepresentation): + def __init__(self, params=None): + super().__init__("SCLSlope") + + def compute_feature(self, scl): + x = _as_float1d(scl) + if x.size < 2: + return np.array(0.0) + t = np.arange(x.size, dtype=np.float64) + return np.array(np.polyfit(t, x, 1)[0]) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class SCLDynamicRange(PhysiologicalRepresentation): + def __init__(self, params=None): + super().__init__("SCLDynamicRange") + + def compute_feature(self, scl): + x = _as_float1d(scl) + if x.size == 0: + return np.array(0.0) + return np.array(np.max(x) - np.min(x)) + + +class _SCRFeature(PhysiologicalRepresentation): + def __init__(self, name, fs=4.0, parameters=None, params=None): + params_dict = parameters or {"fs": [2.0, 4.0, 8.0]} + super().__init__(name, parameters=params_dict, params=params) + self.fs = fs + + def _scr_stats(self, scr): + duration = _as_float1d(scr).size / self.fs + peaks, amplitudes, durations = _detect_scr_peaks(scr, self.fs) + return duration, peaks, amplitudes, durations + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class SCRPeaksPerMinute(_SCRFeature): + def __init__(self, fs=4.0, params=None): + super().__init__("SCRPeaksPerMinute", fs=fs, params=params) + + def compute_feature(self, scr): + duration, peaks, _, _ = self._scr_stats(scr) + if duration <= 0.0: + return np.array(0.0) + return np.array(60.0 * peaks.size / duration) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class SCRAverageAmplitude(_SCRFeature): + def __init__(self, fs=4.0, params=None): + super().__init__("SCRAverageAmplitude", fs=fs, params=params) + + def compute_feature(self, scr): + _, _, amplitudes, _ = self._scr_stats(scr) + if amplitudes.size == 0: + return np.array(0.0) + return np.array(np.mean(amplitudes)) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class SCRAverageDuration(_SCRFeature): + def __init__(self, fs=4.0, params=None): + super().__init__("SCRAverageDuration", fs=fs, params=params) + + def compute_feature(self, scr): + _, _, _, durations = self._scr_stats(scr) + if durations.size == 0: + return np.array(0.0) + return np.array(np.mean(durations)) + + +def _detect_resp_extrema(signal, fs, min_breath_period_s=1.5): + x = _as_float1d(signal) + min_distance = max(1, int(min_breath_period_s * fs)) + if x.size <= min_distance: + empty = np.array([], dtype=int) + return empty, empty + height = np.std(x) * 0.25 + peaks, _ = find_peaks(x, distance=min_distance, height=height) + troughs, _ = find_peaks(-x, distance=min_distance, height=height) + return peaks, troughs + + +def _resp_breath_intervals(signal, fs): + peaks, _ = _detect_resp_extrema(signal, fs) + if peaks.size < 2: + return np.array([], dtype=np.float64) + return np.diff(peaks) / float(fs) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class BreathingRate(PhysiologicalRepresentation): + def __init__(self, fs=500.0, params=None): + super().__init__("BreathingRate") + if params is not None: + fs = params.get("fs", fs) + self.fs = fs + + def compute_feature(self, signal): + intervals = _resp_breath_intervals(signal, self.fs) + if intervals.size == 0: + return np.array(0.0) + return np.array(60.0 / np.mean(intervals)) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class BreathIntervalRMSSD(PhysiologicalRepresentation): + def __init__(self, fs=500.0, params=None): + super().__init__("BreathIntervalRMSSD") + if params is not None: + fs = params.get("fs", fs) + self.fs = fs + + def compute_feature(self, signal): + intervals = _resp_breath_intervals(signal, self.fs) + diffs = _successive_diffs(intervals) + if diffs.size == 0: + return np.array(0.0) + return np.array(np.sqrt(np.mean(diffs**2))) + + +@register_representation([ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class BreathAmplitude(PhysiologicalRepresentation): + def __init__(self, fs=500.0, params=None): + super().__init__("BreathAmplitude") + if params is not None: + fs = params.get("fs", fs) + self.fs = fs + + def compute_feature(self, signal): + x = _as_float1d(signal) + peaks, troughs = _detect_resp_extrema(x, self.fs) + if peaks.size == 0 or troughs.size == 0: + return np.array(0.0) + return np.array(np.mean(x[peaks]) - np.mean(x[troughs])) diff --git a/src/main/python/systemds/scuro/representations/physiological_window.py b/src/main/python/systemds/scuro/representations/physiological_window.py new file mode 100644 index 00000000000..338cfdd2201 --- /dev/null +++ b/src/main/python/systemds/scuro/representations/physiological_window.py @@ -0,0 +1,260 @@ +# ------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +# ------------------------------------------------------------- +import math +import numpy as np + +from systemds.scuro.modality.type import ModalityType +from systemds.scuro.drsearch.operator_registry import register_context_operator +from systemds.scuro.representations.aggregate import Aggregation +from systemds.scuro.representations.representation import RepresentationStats +from systemds.scuro.representations.window_aggregation import ( + Window, + resolve_aggregation_function, + _pad_stack, +) + + +def _estimate_windowed_output_stats( + window_obj, input_stats, estimated_num_windows, representative_window_length +): + in_shape = tuple(int(s) for s in input_stats.output_shape) + if not isinstance(window_obj.aggregation_function, Aggregation): + windowed_input_stats = RepresentationStats( + input_stats.num_instances, (representative_window_length,) + ) + feat_shape = window_obj.aggregation_function.get_output_stats( + windowed_input_stats + ).output_shape + else: + feat_shape = in_shape[1:] + window_obj.stats = RepresentationStats( + input_stats.num_instances, + (estimated_num_windows, *feat_shape), + output_shape_is_known=False, + ) + return window_obj.stats + + +def _estimate_windowed_memory_bytes(window_obj, input_stats): + out_shape = window_obj.get_output_stats(input_stats).output_shape + out_numel = int(np.prod(out_shape)) if len(out_shape) > 0 else 1 + return ( + input_stats.num_instances * out_numel * np.dtype(window_obj.data_type).itemsize + ) + + +def _estimate_windowed_peak_memory_bytes(window_obj, input_stats): + in_shape = tuple(int(s) for s in input_stats.output_shape) + in_numel = int(np.prod(in_shape)) if len(in_shape) > 0 else 1 + input_bytes = ( + input_stats.num_instances * in_numel * np.dtype(window_obj.data_type).itemsize + ) + output_bytes = _estimate_windowed_memory_bytes(window_obj, input_stats) + cpu_peak = int((input_bytes + output_bytes * 2) * 1.2 + 8 * 1024 * 1024) + return {"cpu_peak_bytes": cpu_peak, "gpu_peak_bytes": 0} + + +@register_context_operator([ModalityType.TIMESERIES, ModalityType.PHYSIOLOGICAL]) +class AdaptiveWindow(Window): + granularity_parameter = "base_window_size" + granularity_kind = "length" + + def __init__( + self, + aggregation_function="mean", + base_window_size=256, + overlap=0.5, + min_window_size=64, + params=None, + ): + if params is not None: + aggregation_function = resolve_aggregation_function( + aggregation_function, params + ) + base_window_size = params.get("base_window_size", base_window_size) + overlap = params.get("overlap", overlap) + min_window_size = params.get("min_window_size", min_window_size) + super().__init__("AdaptiveWindow", aggregation_function) + base_window_size = max(1, int(base_window_size)) + min_window_size = max(1, min(int(min_window_size), base_window_size)) + self.parameters.update( + { + "base_window_size": (min(4, base_window_size), base_window_size), + "overlap": [0.25, 0.5, 0.75], + "min_window_size": (min(2, min_window_size), min_window_size), + } + ) + self.base_window_size = base_window_size + self.overlap = overlap + self.min_window_size = min_window_size + + def _estimate_num_windows(self, signal_length): + if signal_length <= 0: + return 1 + step = max(1, int(self.base_window_size * (1 - self.overlap))) + return max(1, math.ceil(signal_length / step)) + + def get_output_stats(self, input_stats: RepresentationStats) -> RepresentationStats: + in_shape = tuple(int(s) for s in input_stats.output_shape) + signal_length = in_shape[0] if in_shape else 0 + return _estimate_windowed_output_stats( + self, + input_stats, + self._estimate_num_windows(signal_length), + self.base_window_size, + ) + + def estimate_output_memory_bytes(self, input_stats: RepresentationStats) -> int: + return _estimate_windowed_memory_bytes(self, input_stats) + + def estimate_peak_memory_bytes(self, input_stats: RepresentationStats) -> dict: + return _estimate_windowed_peak_memory_bytes(self, input_stats) + + def execute(self, modality): + windowed_data = [] + + for signal in modality.data: + local_var = np.array( + [ + np.var(signal[i : i + self.min_window_size]) + for i in range( + 0, len(signal) - self.min_window_size, self.min_window_size + ) + ] + ) + + if local_var.size == 0: + window_sizes = np.array([self.base_window_size]) + else: + norm_var = (local_var - np.min(local_var)) / ( + np.max(local_var) - np.min(local_var) + 1e-6 + ) + window_sizes = np.clip( + self.base_window_size * (1 - 0.5 * norm_var), + self.min_window_size, + self.base_window_size, + ).astype(int) + + windows = [] + start = 0 + while start < len(signal): + current_size = window_sizes[ + min(len(window_sizes) - 1, start // self.min_window_size) + ] + end = min(start + current_size, len(signal)) + window = signal[start:end] + if len(window) > 0: + windows.append(window) + start += max(1, int(current_size * (1 - self.overlap))) + + processed_windows = _pad_stack( + [self.aggregation_function.compute_feature(w) for w in windows] + ) + windowed_data.append(processed_windows) + + return _pad_stack(windowed_data) + + +@register_context_operator([ModalityType.TIMESERIES, ModalityType.PHYSIOLOGICAL]) +class PhysiologicalEventWindow(Window): + granularity_parameter = "min_distance" + granularity_kind = "length" + + def __init__( + self, + aggregation_function="mean", + event_threshold=0.5, + min_distance=100, + params=None, + ): + if params is not None: + aggregation_function = resolve_aggregation_function( + aggregation_function, params + ) + event_threshold = params.get("event_threshold", event_threshold) + min_distance = params.get("min_distance", min_distance) + super().__init__("PhysiologicalEventWindow", aggregation_function) + min_distance = max(1, int(min_distance)) + self.parameters.update( + { + "event_threshold": [0.3, 0.5, 0.7], + "min_distance": (min(4, min_distance), min_distance), + } + ) + self.event_threshold = event_threshold + self.min_distance = min_distance + + def _estimate_num_windows(self, signal_length): + if signal_length <= 0: + return 1 + return max(1, math.ceil(signal_length / max(1, self.min_distance))) + + def get_output_stats(self, input_stats: RepresentationStats) -> RepresentationStats: + in_shape = tuple(int(s) for s in input_stats.output_shape) + signal_length = in_shape[0] if in_shape else 0 + return _estimate_windowed_output_stats( + self, + input_stats, + self._estimate_num_windows(signal_length), + self.min_distance, + ) + + def estimate_output_memory_bytes(self, input_stats: RepresentationStats) -> int: + return _estimate_windowed_memory_bytes(self, input_stats) + + def estimate_peak_memory_bytes(self, input_stats: RepresentationStats) -> dict: + return _estimate_windowed_peak_memory_bytes(self, input_stats) + + def execute(self, modality): + windowed_data = [] + + for signal in modality.data: + normalized = (signal - np.mean(signal)) / (np.std(signal) + 1e-6) + + peaks = [] + last_peak = -self.min_distance + for i in range(1, len(normalized) - 1): + if ( + normalized[i] > self.event_threshold + and normalized[i] > normalized[i - 1] + and normalized[i] > normalized[i + 1] + and i - last_peak >= self.min_distance + ): + peaks.append(i) + last_peak = i + + windows = [signal[peaks[i] : peaks[i + 1]] for i in range(len(peaks) - 1)] + + if not windows: + windows = [ + w + for w in np.array_split( + signal, max(1, len(signal) // self.min_distance) + ) + if len(w) > 0 + ] + + processed_windows = _pad_stack( + [self.aggregation_function.compute_feature(w) for w in windows] + ) + windowed_data.append(processed_windows) + + return _pad_stack(windowed_data) diff --git a/src/main/python/systemds/scuro/representations/sum.py b/src/main/python/systemds/scuro/representations/sum.py index 4f658020f1e..58cac2e7210 100644 --- a/src/main/python/systemds/scuro/representations/sum.py +++ b/src/main/python/systemds/scuro/representations/sum.py @@ -41,7 +41,7 @@ def __init__(self, params=None): self.needs_alignment = True def execute(self, modalities: List[Modality]): - data = np.asarray( + data = np.array( modalities[0].data, dtype=modalities[0].metadata[0]["data_layout"]["type"], ) @@ -54,31 +54,28 @@ def execute(self, modalities: List[Modality]): return data def get_output_stats(self, input_stats_list) -> RepresentationStats: - if isinstance(input_stats_list, RepresentationStats): - return input_stats_list - - stats_list = list(input_stats_list) + stats_list = self._fusion_input_stats(input_stats_list) if not stats_list: return RepresentationStats(0, (0,)) - max_dim = max([stats.output_shape[-1] for stats in stats_list]) - return RepresentationStats(stats_list[0].num_instances, (max_dim,)) - - def estimate_peak_memory_bytes(self, input_stats_list) -> dict: - elem_size = np.dtype(np.float64).itemsize - - def stats_payload_bytes(s: RepresentationStats) -> int: - numel = int(np.prod(s.output_shape)) if len(s.output_shape) > 0 else 1 - return int(s.num_instances * numel * elem_size) + num_instances = max(s.num_instances for s in stats_list) + rank = len(stats_list[0].output_shape) + if rank > 0 and all(len(s.output_shape) == rank for s in stats_list): + output_shape = tuple( + max(s.output_shape[d] for s in stats_list) for d in range(rank) + ) + else: + output_shape = max(stats_list, key=self._stats_num_elements).output_shape + output_shape_is_known = all(s.output_shape_is_known for s in stats_list) + return RepresentationStats(num_instances, output_shape, output_shape_is_known) - first_bytes = stats_payload_bytes(input_stats_list[0]) - max_other_bytes = 0 - if len(input_stats_list) > 1: - max_other_bytes = max(stats_payload_bytes(s) for s in input_stats_list[1:]) + def estimate_peak_memory_bytes(self, input_stats) -> dict: + stats_list = self._as_stats_list(input_stats) + input_bytes = sum(self._stats_bytes(s) for s in stats_list) + output_bytes = self._stats_bytes(self.get_output_stats(input_stats)) - ufunc_workspace_bytes = int(0.1 * max(first_bytes, max_other_bytes)) - cpu_peak = int( - (first_bytes + max_other_bytes + ufunc_workspace_bytes) * 1.15 - + 8 * 1024 * 1024 + raw_bytes = self._raw_input_bytes(input_stats) + cpu_peak = ( + int((raw_bytes + input_bytes + output_bytes) * 1.15) + 8 * 1024 * 1024 ) return {"cpu_peak_bytes": cpu_peak, "gpu_peak_bytes": 0} diff --git a/src/main/python/systemds/scuro/representations/tabular_features.py b/src/main/python/systemds/scuro/representations/tabular_features.py new file mode 100644 index 00000000000..9734faa39db --- /dev/null +++ b/src/main/python/systemds/scuro/representations/tabular_features.py @@ -0,0 +1,59 @@ +# ------------------------------------------------------------- +# +# Licensed to the Apache Software Foundation (ASF) under one +# or more contributor license agreements. See the NOTICE file +# distributed with this work for additional information +# regarding copyright ownership. The ASF licenses this file +# to you under the Apache License, Version 2.0 (the +# "License"); you may not use this file except in compliance +# with the License. You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, +# software distributed under the License is distributed on an +# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +# KIND, either express or implied. See the License for the +# specific language governing permissions and limitations +# under the License. +# +# ------------------------------------------------------------- +import numpy as np + +from systemds.scuro.dataloader.tabular_loader import TabularStats +from systemds.scuro.modality.transformed import TransformedModality +from systemds.scuro.modality.type import ModalityType +from systemds.scuro.drsearch.operator_registry import register_representation +from systemds.scuro.representations.representation import RepresentationStats +from systemds.scuro.representations.unimodal import UnimodalRepresentation + + +@register_representation(ModalityType.EMBEDDING) +class TabularFeatures(UnimodalRepresentation): + def __init__(self, params=None): + super().__init__("TabularFeatures", ModalityType.EMBEDDING, None) + self.data_type = np.float32 + + def get_output_stats(self, input_stats: TabularStats) -> RepresentationStats: + return RepresentationStats( + input_stats.num_instances, input_stats.output_shape, dtype=self.data_type + ) + + def estimate_output_memory_bytes(self, input_stats: TabularStats) -> int: + return ( + input_stats.num_instances + * input_stats.num_features + * np.dtype(self.data_type).itemsize + ) + + def estimate_peak_memory_bytes(self, input_stats: TabularStats) -> dict: + return { + "cpu_peak_bytes": self.estimate_output_memory_bytes(input_stats) * 2, + "gpu_peak_bytes": 0, + } + + def transform(self, modality, params=None): + transformed_modality = TransformedModality(modality, self) + transformed_modality.data_type = self.data_type + transformed_modality.data = np.array(modality.data, dtype=self.data_type) + return transformed_modality diff --git a/src/main/python/systemds/scuro/representations/text_context_with_indices.py b/src/main/python/systemds/scuro/representations/text_context_with_indices.py index 4de53698d7a..6d8964c1d14 100644 --- a/src/main/python/systemds/scuro/representations/text_context_with_indices.py +++ b/src/main/python/systemds/scuro/representations/text_context_with_indices.py @@ -177,6 +177,7 @@ def get_output_stats(self, input_stats: TextStats) -> RepresentationStats: ), self.max_words, ), + dtype=self.data_type, ) def estimate_output_memory_bytes(self, input_stats: TextStats) -> int: @@ -318,6 +319,7 @@ def get_output_stats(self, input_stats: TextStats) -> RepresentationStats: ), self.max_words, ), + dtype=self.data_type, ) def estimate_output_memory_bytes(self, input_stats: TextStats) -> int: diff --git a/src/main/python/systemds/scuro/representations/tfidf.py b/src/main/python/systemds/scuro/representations/tfidf.py index 3c3d894c173..bea18a56024 100644 --- a/src/main/python/systemds/scuro/representations/tfidf.py +++ b/src/main/python/systemds/scuro/representations/tfidf.py @@ -47,7 +47,10 @@ def get_output_stats(self, input_stats: TextStats) -> RepresentationStats: 100_000, max(1000, input_stats.num_instances * input_stats.max_length) ) return RepresentationStats( - input_stats.num_instances, (vocab_estimate,), output_shape_is_known=False + input_stats.num_instances, + (vocab_estimate,), + output_shape_is_known=False, + dtype=self.data_type, ) def estimate_output_memory_bytes(self, input_stats: TextStats) -> int: diff --git a/src/main/python/systemds/scuro/representations/timeseries_representations.py b/src/main/python/systemds/scuro/representations/timeseries_representations.py index 6bf7f38d132..2c7fbbe7404 100644 --- a/src/main/python/systemds/scuro/representations/timeseries_representations.py +++ b/src/main/python/systemds/scuro/representations/timeseries_representations.py @@ -23,28 +23,75 @@ from systemds.scuro.modality.type import ModalityType from systemds.scuro.modality.transformed import TransformedModality -from systemds.scuro.representations.representation import RepresentationStats +from systemds.scuro.representations.representation import ( + CONTAINER_ARRAY, + CONTAINER_LIST, + RepresentationStats, +) from systemds.scuro.representations.unimodal import UnimodalRepresentation +from systemds.scuro.representations.utils import dense_instance_batch from systemds.scuro.drsearch.operator_registry import ( register_representation, register_context_representation_operator, ) +import warnings + +warnings.filterwarnings( + "ignore", + message=r"Precision loss occurred in moment calculation", + category=RuntimeWarning, +) + class TimeSeriesRepresentation(UnimodalRepresentation): - def __init__(self, name, parameters=None, params=None): + + def __init__( + self, + name, + parameters=None, + params=None, + self_contained=False, + min_input_length=1, + ): if params is None: params = {} - - super().__init__(name, ModalityType.EMBEDDING, parameters, False) + self.min_input_length = min_input_length + super().__init__(name, ModalityType.EMBEDDING, parameters, self_contained) + + @staticmethod + def _input_length(input_stats) -> int: + shape = getattr(input_stats, "output_shape", ()) or () + return int(shape[0]) if shape else 0 + + def check_preconditions(self, input_stats): + length = self._input_length(input_stats) + if length < self.min_input_length: + return ( + f"{self.name} needs >= {self.min_input_length} samples, " + f"input provides {length}" + ) + return None def compute_feature(self, signal): raise NotImplementedError("Subclasses should implement this method.") + def compute_features_batched(self, data): + return np.asarray(self.compute_feature(data, axis=-1)) + def transform(self, modality, aggregation=None): transformed_modality = TransformedModality( modality, self, self.output_modality_type ) + dtype = modality.metadata[0]["data_layout"]["type"] + batch = dense_instance_batch(modality.data) + if batch is not None: + features = self.compute_features_batched(batch) + if features.ndim == 1: + features = features[:, None] + transformed_modality.data = features.astype(dtype) + return transformed_modality + result = [] for signal in modality.data: @@ -56,23 +103,45 @@ def transform(self, modality, aggregation=None): np.pad(r, (0, maxlen - r.size), mode="constant", constant_values=0.0) for r in result ] - dtype = modality.metadata[0]["data_layout"]["type"] transformed_modality.data = np.vstack(np.asarray(padded_result)).astype(dtype) return transformed_modality def get_output_stats(self, input_stats): - return RepresentationStats(input_stats.num_instances, (1,)) + return RepresentationStats( + input_stats.num_instances, (1,), input_stats.output_shape_is_known + ) + + @staticmethod + def _num_elements(shape) -> int: + n = 1 + for d in shape: + n *= int(d) + return n def estimate_output_memory_bytes(self, input_stats): - # TODO: adapt this to the actual output shapes and transformations - return input_stats.num_instances * 4 + out_stats = self.get_output_stats(input_stats) + return ( + int(out_stats.num_instances) + * self._num_elements(out_stats.output_shape) + * np.dtype(np.float32).itemsize + ) def estimate_peak_memory_bytes(self, input_stats): - # TODO: adapt this to the actual output shapes and transformations - return { - "cpu_peak_bytes": self.estimate_output_memory_bytes(input_stats), - "gpu_peak_bytes": 0, - } + input_bytes = ( + int(input_stats.num_instances) + * self._num_elements(input_stats.output_shape) + * np.dtype(np.float32).itemsize + ) + output_bytes = self.estimate_output_memory_bytes(input_stats) + batch_bytes = ( + input_bytes + if getattr(input_stats, "container", CONTAINER_ARRAY) == CONTAINER_LIST + else 0 + ) + cpu_peak = ( + int((input_bytes + batch_bytes + 3 * output_bytes) * 1.15) + 4 * 1024 * 1024 + ) + return {"cpu_peak_bytes": cpu_peak, "gpu_peak_bytes": 0} @register_representation([ModalityType.TIMESERIES, ModalityType.PHYSIOLOGICAL]) @@ -126,7 +195,7 @@ def compute_feature(self, signal, axis=-1): @register_context_representation_operator(ModalityType.AUDIO) class Std(TimeSeriesRepresentation): def __init__(self, params=None): - super().__init__("Std") + super().__init__("Std", min_input_length=2) def compute_feature(self, signal, axis=-1): return np.array(np.std(signal, axis=axis)) @@ -138,7 +207,7 @@ def compute_feature(self, signal, axis=-1): @register_context_representation_operator(ModalityType.AUDIO) class Skew(TimeSeriesRepresentation): def __init__(self, params=None): - super().__init__("Skew") + super().__init__("Skew", min_input_length=3) def compute_feature(self, signal, axis=-1): return np.array(stats.skew(signal, axis=axis)) @@ -152,17 +221,36 @@ def __init__(self, quantile=0.9, params=None): super().__init__( "Qunatile", {"quantile": [0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99]} ) + if params is not None: + quantile = params.get("quantile", quantile) self.quantile = quantile def compute_feature(self, signal, axis=-1): return np.array(np.quantile(signal, self.quantile, axis=axis)) + def compute_features_batched(self, data): + features = np.asarray(np.quantile(data, self.quantile, axis=-1)) + if np.ndim(self.quantile) == 0: + return features + + return np.moveaxis(features, 0, -1) + + def get_output_stats(self, input_stats): + n_quantiles = np.atleast_1d(self.quantile).size + return RepresentationStats( + input_stats.num_instances, + (n_quantiles,), + input_stats.output_shape_is_known, + ) + @register_representation([ModalityType.TIMESERIES, ModalityType.PHYSIOLOGICAL]) @register_context_representation_operator(ModalityType.TIMESERIES) @register_context_representation_operator(ModalityType.PHYSIOLOGICAL) @register_context_representation_operator(ModalityType.AUDIO) class Kurtosis(TimeSeriesRepresentation): + min_input_length = 4 # the fourth moment is undefined below four samples + def __init__(self, params=None): super().__init__("Kurtosis") @@ -187,22 +275,89 @@ def compute_feature(self, signal, axis=-1): @register_context_representation_operator(ModalityType.PHYSIOLOGICAL) class ZeroCrossingRate(TimeSeriesRepresentation): def __init__(self, params=None): - super().__init__("ZeroCrossingRate") + super().__init__("ZeroCrossingRate", min_input_length=2) def compute_feature(self, signal, axis=-1): return np.array(np.sum(np.diff(np.signbit(signal), axis=axis) != 0, axis=axis)) +@register_representation([ModalityType.TIMESERIES, ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.TIMESERIES) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class LastValue(TimeSeriesRepresentation): + def __init__(self, params=None): + super().__init__("LastValue") + + def compute_feature(self, signal, axis=-1): + return np.take(signal, -1, axis=axis) + + +@register_representation([ModalityType.TIMESERIES, ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.TIMESERIES) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class TransitionCount(TimeSeriesRepresentation): + def __init__(self, threshold=0.0, params=None): + super().__init__( + "TransitionCount", + parameters={"threshold": [0.0, 1.0, 5.0, 10.0]}, + min_input_length=2, + ) + if params is not None: + threshold = params.get("threshold", threshold) + self.threshold = threshold + + def compute_feature(self, signal, axis=-1): + return np.array( + np.sum(np.abs(np.diff(signal, axis=axis)) > self.threshold, axis=axis) + ) + + +@register_representation([ModalityType.TIMESERIES, ModalityType.PHYSIOLOGICAL]) +@register_context_representation_operator(ModalityType.TIMESERIES) +@register_context_representation_operator(ModalityType.PHYSIOLOGICAL) +class ObservationDensity(TimeSeriesRepresentation): + def __init__(self, params=None): + super().__init__("ObservationDensity", min_input_length=2) + + def compute_feature(self, signal, axis=-1): + n = signal.shape[axis] + transitions = np.sum(np.diff(signal, axis=axis) != 0, axis=axis) + return np.array((transitions + 1) / n) + + @register_representation([ModalityType.TIMESERIES, ModalityType.PHYSIOLOGICAL]) @register_context_representation_operator(ModalityType.TIMESERIES) @register_context_representation_operator(ModalityType.PHYSIOLOGICAL) class ACF(TimeSeriesRepresentation): def __init__(self, k=1, params=None): - super().__init__("ACF", {"k": [1, 2, 5, 10, 20, 25, 50, 100, 200, 500]}) + super().__init__( + "ACF", {"k": [1, 2, 5, 10, 20, 25, 50, 100, 200, 500]}, min_input_length=2 + ) if params is not None: k = params.get("k", k) self.k = k + def filter_parameter_domain(self, name, values, input_stats): + if name != "k": + return values + length = self._input_length(input_stats) + if length <= 1: + return values + usable = [k for k in values if 0 < int(k) < length] + return usable or [1] + + def check_preconditions(self, input_stats): + failure = super().check_preconditions(input_stats) + if failure: + return failure + length = self._input_length(input_stats) + if int(self.k) >= length: + return ( + f"ACF lag k={int(self.k)} needs > {int(self.k)} samples, " + f"input provides {length}" + ) + return None + def compute_feature(self, signal, axis=-1): x = np.asarray(signal, dtype=np.float64) x = x - np.mean(x, axis=axis, keepdims=True) @@ -235,22 +390,39 @@ def get_k_values(self, max_length, percent=0.2, num=10, log=False): @register_context_representation_operator(ModalityType.TIMESERIES) @register_context_representation_operator(ModalityType.PHYSIOLOGICAL) class FrequencyMagnitude(TimeSeriesRepresentation): - def __init__(self, params=None): - super().__init__("FrequencyMagnitude") + def __init__(self, params=None, self_contained=True): + super().__init__("FrequencyMagnitude", min_input_length=2) def compute_feature(self, signal, axis=-1): return np.array(np.abs(np.fft.rfft(signal, axis=axis))) + def get_output_stats(self, input_stats): + n = self._num_elements(input_stats.output_shape) + out_len = n // 2 + 1 if n > 0 else 0 + return RepresentationStats( + input_stats.num_instances, (out_len,), input_stats.output_shape_is_known + ) + @register_representation([ModalityType.TIMESERIES, ModalityType.PHYSIOLOGICAL]) @register_context_representation_operator(ModalityType.PHYSIOLOGICAL) @register_context_representation_operator(ModalityType.TIMESERIES) class SpectralCentroid(TimeSeriesRepresentation): def __init__(self, fs=1.0, params=None): - super().__init__("SpectralCentroid", parameters={"fs": [0.5, 1.0, 2.0]}) + super().__init__("SpectralCentroid", min_input_length=2) if params is not None: fs = params.get("fs", fs) - self.fs = fs + self.fs = float(fs) + + def get_current_parameters(self): + current_params = super().get_current_parameters() + current_params["fs"] = self.fs + return current_params + + def configure_for_input(self, input_stats): + sampling_rate = getattr(input_stats, "sampling_rate", None) + if sampling_rate: + self.fs = float(sampling_rate) def compute_feature(self, signal, axis=-1): signal = np.asarray(signal, dtype=np.float64) @@ -270,22 +442,42 @@ def compute_feature(self, signal, axis=-1): @register_context_representation_operator(ModalityType.PHYSIOLOGICAL) @register_context_representation_operator(ModalityType.TIMESERIES) class BandpowerFFT(TimeSeriesRepresentation): - def __init__(self, fs=1.0, f1=0.0, f2=0.5, params=None): + def __init__(self, fs=1.0, band_low=0.0, band_width=0.5, params=None): super().__init__( "BandpowerFFT", - parameters={"fs": [0.5, 1.0], "f1": [0.0, 1.0], "f2": [0.5, 1.0]}, + parameters={ + "band_low": [0.0, 0.25, 0.5], + "band_width": [0.25, 0.5, 1.0], + }, + min_input_length=2, ) if params is not None: fs = params.get("fs", fs) - f1 = params.get("f1", f1) - f2 = params.get("f2", f2) - self.fs = fs - self.f1 = f1 - self.f2 = f2 + band_low = params.get("band_low", band_low) + band_width = params.get("band_width", band_width) + self.fs = float(fs) + self.band_low = float(band_low) + self.band_width = float(band_width) + + @property + def band_high(self) -> float: + return min(1.0, self.band_low + self.band_width) + + def get_current_parameters(self): + current_params = super().get_current_parameters() + current_params["fs"] = self.fs # bound to the data, not searched + return current_params + + def configure_for_input(self, input_stats): + sampling_rate = getattr(input_stats, "sampling_rate", None) + if sampling_rate: + self.fs = float(sampling_rate) def compute_feature(self, signal, axis=-1): signal = np.asarray(signal, dtype=np.float64) n = signal.shape[axis] + nyquist = self.fs / 2.0 + self.f1, self.f2 = self.band_low * nyquist, self.band_high * nyquist frequency_magnitude = FrequencyMagnitude().compute_feature(signal, axis=axis) frequencies = np.fft.rfftfreq(n, d=1.0 / self.fs) diff --git a/src/main/python/systemds/scuro/representations/utils.py b/src/main/python/systemds/scuro/representations/utils.py index 7551c6cb2bf..5041e18770c 100644 --- a/src/main/python/systemds/scuro/representations/utils.py +++ b/src/main/python/systemds/scuro/representations/utils.py @@ -24,6 +24,29 @@ import numpy as np +def dense_instance_batch(data): + if isinstance(data, np.ndarray): + if data.ndim == 2 and np.issubdtype(data.dtype, np.number): + return data + return None + + if not isinstance(data, (list, tuple)) or len(data) == 0: + return None + + first = data[0] + if ( + not isinstance(first, np.ndarray) + or first.ndim != 1 + or not np.issubdtype(first.dtype, np.number) + ): + return None + for instance in data: + if not isinstance(instance, np.ndarray) or instance.shape != first.shape: + return None + + return np.asarray(data) + + def pad_sequences(sequences, maxlen=None, dtype="float32", value=0): if maxlen is None: maxlen = max([len(seq) for seq in sequences]) diff --git a/src/main/python/systemds/scuro/representations/wav2vec.py b/src/main/python/systemds/scuro/representations/wav2vec.py index 5e03baf8bc4..c9f4025579a 100644 --- a/src/main/python/systemds/scuro/representations/wav2vec.py +++ b/src/main/python/systemds/scuro/representations/wav2vec.py @@ -37,14 +37,34 @@ @register_representation(ModalityType.AUDIO) class Wav2Vec(UnimodalRepresentation): + cache_in_worker = True + instance_parallel = True + + MODEL_NAME = "facebook/wav2vec2-base-960h" + def __init__(self, params=None): super().__init__("Wav2Vec", ModalityType.TIMESERIES, {}) - self.processor = Wav2Vec2Processor.from_pretrained( - "facebook/wav2vec2-base-960h" - ) - self.model = Wav2Vec2Model.from_pretrained( - "facebook/wav2vec2-base-960h" - ).float() + self._processor = None + self._model = None + + @staticmethod + def _from_pretrained(loader_cls, name): + try: + return loader_cls.from_pretrained(name, local_files_only=True) + except Exception: + return loader_cls.from_pretrained(name) + + @property + def processor(self): + if self._processor is None: + self._processor = self._from_pretrained(Wav2Vec2Processor, self.MODEL_NAME) + return self._processor + + @property + def model(self): + if self._model is None: + self._model = self._from_pretrained(Wav2Vec2Model, self.MODEL_NAME).float() + return self._model def transform(self, modality, aggregation=None): transformed_modality = TransformedModality( diff --git a/src/main/python/systemds/scuro/representations/window_aggregation.py b/src/main/python/systemds/scuro/representations/window_aggregation.py index a9a1f1eb41b..713a398312a 100644 --- a/src/main/python/systemds/scuro/representations/window_aggregation.py +++ b/src/main/python/systemds/scuro/representations/window_aggregation.py @@ -29,11 +29,25 @@ from systemds.scuro.representations.aggregate import Aggregation from systemds.scuro.representations.context import Context from systemds.scuro.representations.representation import ( + CONTAINER_ARRAY, + CONTAINER_LIST, NDARRAY_OBJECT_OVERHEAD_BYTES, Representation, RepresentationStats, stats_itemsize, ) +from systemds.scuro.representations.utils import dense_instance_batch + +_ACCEPTS_AXIS_CACHE = {} + + +def _accepts_axis(compute_feature): + func = getattr(compute_feature, "__func__", compute_feature) + accepts = _ACCEPTS_AXIS_CACHE.get(func) + if accepts is None: + accepts = "axis" in inspect.signature(compute_feature).parameters + _ACCEPTS_AXIS_CACHE[func] = accepts + return accepts def nested_aggregation_param_names(agg_cls): @@ -97,6 +111,23 @@ def _append_tail_row(full_result, tail_result): return np.concatenate([full_result, tail_row[None, ...]]) +def _append_tail_rows(full_result, tail_result): + full_result = np.asarray(full_result) + tail_result = np.asarray(tail_result) + target_shape = full_result.shape[2:] + if tail_result.shape[1:] == target_shape: + tail_rows = tail_result + else: + tail_rows = np.zeros( + (full_result.shape[0], *target_shape), dtype=full_result.dtype + ) + slices = tuple( + slice(0, min(d, s)) for d, s in zip(target_shape, tail_result.shape[1:]) + ) + tail_rows[(slice(None), *slices)] = tail_result[(slice(None), *slices)] + return np.concatenate([full_result, tail_rows[:, None, ...]], axis=1) + + def resolve_aggregation_function(aggregation_function, params): if params is None: return aggregation_function @@ -300,10 +331,32 @@ def estimate_peak_memory_bytes(self, input_stats: RepresentationStats) -> dict: + padded_elems * stats_itemsize(input_stats) ) - cpu_peak = int((input_bytes + list_bytes + pad_bytes) * 1.15 + 16 * 1024 * 1024) + batch_bytes = ( + input_bytes + if getattr(input_stats, "container", CONTAINER_ARRAY) == CONTAINER_LIST + else 0 + ) + + cpu_peak = int( + (input_bytes + batch_bytes + list_bytes + pad_bytes) * 1.15 + + 16 * 1024 * 1024 + ) return {"cpu_peak_bytes": cpu_peak, "gpu_peak_bytes": 0} def execute(self, modality): + batch = self._dense_batch(modality) + if batch is not None: + windowed_data = self.window_aggregate_single_level_batched(batch) + if windowed_data is not None: + if self.pad: + data_type = modality.metadata[0]["data_layout"]["type"] + if data_type != "str": + windowed_data = windowed_data.astype(data_type) + else: + windowed_data = list(windowed_data) + self.assert_output_stats(windowed_data) + return windowed_data + windowed_data = [] original_lengths = [] for instance in modality.data: @@ -352,6 +405,56 @@ def execute(self, modality): self.assert_output_stats(windowed_data) return windowed_data + def _dense_batch(self, modality): + if modality.get_data_layout() != DataLayout.SINGLE_LEVEL: + return None + if not _accepts_axis(self.aggregation_function.compute_feature): + return None + + batch = dense_instance_batch(modality.data) + if batch is None: + return None + + batch = batch.view() + batch.setflags(write=False) + return batch + + def window_aggregate_single_level_batched(self, data): + num_instances, length = data.shape + new_length = math.ceil(length / self.window_size) + cut_length = (new_length - 1) * self.window_size + tail = data[:, cut_length:] + compute_feature = self.aggregation_function.compute_feature + + if new_length <= 1: + if not tail.shape[1]: + raise ValueError( + "Cannot window-aggregate an empty instance " + f"(window_size={self.window_size})." + ) + if tail.shape[1] < self.window_size: + pad_len = self.window_size - tail.shape[1] + tail = np.pad(tail, ((0, 0), (0, pad_len)), mode="constant") + result = np.asarray(compute_feature(tail, axis=1)) + if result.shape[0] != num_instances: + return None + return result[:, None, ...] + + full_batches = data[:, :cut_length].reshape( + num_instances, new_length - 1, self.window_size + ) + full_result = np.asarray(compute_feature(full_batches, axis=2)) + if full_result.shape[:2] != (num_instances, new_length - 1): + return None + + if tail.shape[1]: + tail_result = np.asarray(compute_feature(tail, axis=1)) + if tail_result.shape[0] != num_instances: + return None + full_result = _append_tail_rows(full_result, tail_result) + + return full_result + def window_aggregate_single_level(self, instance, new_length): if isinstance(instance, str): return instance @@ -359,7 +462,7 @@ def window_aggregate_single_level(self, instance, new_length): arr = np.asarray(instance) cut_length = (new_length - 1) * self.window_size tail = arr[cut_length:] - sig = inspect.signature(self.aggregation_function.compute_feature) + takes_axis = _accepts_axis(self.aggregation_function.compute_feature) if new_length <= 1: if not tail.size: raise ValueError( @@ -374,7 +477,7 @@ def window_aggregate_single_level(self, instance, new_length): pad_width = [(0, 0)] * tail.ndim pad_width[0] = (0, pad_len) tail = np.pad(tail, pad_width=pad_width, mode="constant") - if "axis" in sig.parameters: + if takes_axis: return np.array([self.aggregation_function.compute_feature(tail)]) tail_result = self.aggregation_function.compute_feature(tail) return ( @@ -386,7 +489,7 @@ def window_aggregate_single_level(self, instance, new_length): new_length - 1, self.window_size, *arr.shape[1:] ) - if "axis" in sig.parameters: + if takes_axis: full_result = self.aggregation_function.compute_feature( full_batches, axis=1 ) @@ -495,8 +598,7 @@ def execute(self, modality): self.num_windows, window_size, *instance.shape[1:] ) - sig = inspect.signature(self.aggregation_function.compute_feature) - if "axis" in sig.parameters: + if _accepts_axis(self.aggregation_function.compute_feature): f = self.aggregation_function.compute_feature(full_batches, axis=1) else: f = np.stack( diff --git a/src/main/python/systemds/scuro/representations/word2vec.py b/src/main/python/systemds/scuro/representations/word2vec.py index bc1c8791f20..fd1e148e117 100644 --- a/src/main/python/systemds/scuro/representations/word2vec.py +++ b/src/main/python/systemds/scuro/representations/word2vec.py @@ -58,7 +58,9 @@ def __init__(self, vector_size=150, min_count=1, output_file=None, params=None): self.data_type = np.float32 def get_output_stats(self, input_stats: TextStats) -> RepresentationStats: - return RepresentationStats(input_stats.num_instances, (self.vector_size,)) + return RepresentationStats( + input_stats.num_instances, (self.vector_size,), dtype=self.data_type + ) def estimate_output_memory_bytes(self, input_stats: TextStats) -> int: return ( diff --git a/src/main/python/tests/scuro/test_operator_registry.py b/src/main/python/tests/scuro/test_operator_registry.py index 443cc039d6b..93afba342b0 100644 --- a/src/main/python/tests/scuro/test_operator_registry.py +++ b/src/main/python/tests/scuro/test_operator_registry.py @@ -62,6 +62,9 @@ Quantile, ZeroCrossingRate, FrequencyMagnitude, + LastValue, + TransitionCount, + ObservationDensity, ) from systemds.scuro.modality.type import ModalityType from systemds.scuro.representations.average import Average @@ -74,6 +77,10 @@ from systemds.scuro.representations.hadamard import Hadamard from systemds.scuro.representations.resnet import ResNet from systemds.scuro.representations.multimodal_attention_fusion import AttentionFusion +from systemds.scuro.representations.physiological_window import ( + AdaptiveWindow, + PhysiologicalEventWindow, +) class TestOperatorRegistry(unittest.TestCase): @@ -113,6 +120,9 @@ def test_timeseries_representations_in_registry(self): Kurtosis, RMS, ZeroCrossingRate, + LastValue, + TransitionCount, + ObservationDensity, ACF, FrequencyMagnitude, SpectralCentroid, @@ -132,6 +142,8 @@ def test_context_operator_in_registry(self): WindowAggregation, StaticWindow, DynamicWindow, + AdaptiveWindow, + PhysiologicalEventWindow, ] assert registry.get_context_operators(ModalityType.TEXT) == [ SentenceBoundarySplitIndices,