The resample class introduced in #38 has some support for resampling index values together with the data.

But, we still don't handle some cases:

• We don't check if the index in actually increasing
• Handling of datetime index might be possible (added in 82a967c)

We need to decide what the input for each algorithm should be.

The entire pipeline needs to be capable of processing full day recordings. However, realistically, only the GSD algorithm needs to see the full data.

All algorithms afterwards only process individual GSs. Hence, I think it makes sense to have just the IMU data of 1 GS as the minimal input of these algorithm. But, this means we need to write a wrapper that makes it easy to run an algorithm on multiple GSs.

We need to define how the final result data structures should look like.

In general I would be in favor of a simple pd.DataFrame, however, while the summary parameters are just a single value per WB, we also need to store the stride level parameter, the per-second parameter, the turn-level parameter and potentially more information is strange data structure.

Once decided we should also update our data loader for reference data to prduce this ouput formar

At the moment we only test 3.10 on Linux to save on CI minutes. Once development has slowed down, we should test all supported Python versions and operating systems

Many algorithms require to resample/downsample the data to 40 Hz.

The original implementation uses "simple" linear interpolation to do that. This might introduce artifacts, as the interpolation is done without prior lowpass filtering.

As this shows, resampling is a difficult topic that we should solve with centralized utilities.

For this we should:

• research different approaches already implemented in scipy (resample vs. resample_poly vs. decimate)
• Test all approaches on some example data and make a decision, if we see any clear advantages or disadvantages of any of these approaches.
• Implement a resampling as part of the data_transform module

This issue is meant to track the implementation of the {CADB_CADC} in python

Short description: {The Cadence Detection toolbox contains code (MATLAB, R2018b) for estimating the cadence during each detected gait sequence. The method is based on steps detection and demarcation. Step duration is defined as the period between two consecutive feet initial contacts (consecutive heel strikes from left and right foot). Then, the instantaneous cadence is estimated as the inverse function of the step duration (in minute unit to have steps/min unit for the cadence). Eventually, for each algorithm, the mean value of the cadence over each gait sequence is computed. The algorithm provides the instantaneous cadence at every second during the gait sequence.

For detection of steps two algorithms are available. Both use the tri-axial acceleration norm after pre-processing to filter the noise and enhance the step-related pattern.

Pre-Processing includes: detrending, low-pass filtering (FIR, fc=3.2 Hz), smoothing using Savitsky-Golay filter (order=7, frame length =21) and Gaussian smoothing, and continuous wavelet transform (CWT, scale 10, ‘gaus2’) to enhance the relevant steps-related features (peaks) in acceleration signal.

Steps detection and demarcation (and subsequently the cadence) can be performed by selecting one of the two implemented algorithms:

With the algorithm ‘HKLee_imp’ (implemented by function HKLee_improved.m), the pre-processed signal is further processed using morphological filters, according to methods described in [1,2]. Finally, the timing of steps, used to estimate step frequency, are detected as maxima between successive zero-crossings.

With the algorithm ‘Shin_imp’ (implemented by function Shin_improved.m), the timing of steps, used to estimate step frequency, are identified as the zero-crossing on the positive slopes of pre-processed signal [1, 3].

The cadence detection algorithm is implemented in the main function CADENCE.m. The script example_GSD_CAD.m contains an exemplary application of the algorithm, using input data in the standardized format adopted by MobiliseD project [4]. The algorithm requires also the start and end of each gait sequence, in the format provided by the Gait Sequence detection algorithm.

Note that in the validation study [3] cadence estimation using ‘HKLee_imp’ is referred as CADB, and using ‘Shin_imp’ as CADc.}

Link to original implementation: {Link}

Reference Papers: {[1] Paraschiv-Ionescu, A, Soltani A, and Aminian K. "Real-world speed estimation using single trunk IMU: methodological challenges for impaired gait patterns." 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). IEEE, 2020.

[2] Lee, H-K., et al. "Computational methods to detect step events for normal and pathological gait evaluation using accelerometer." Electronics letters 46.17 (2010): 1185-1187.

[3] [1] Shin, Seung Hyuck, and Chan Gook Park. "Adaptive step length estimation algorithm using optimal parameters and movement status awareness”. Medical engineering & physics 33.9 (2011): 1064-1071.

[4] Micó-Amigo, M. E., Bonci, T., Paraschiv-Ionescu, A., Ullrich, M., Kirk, C., Soltani, A., ... & Del Din, S. (2022). Assessing real-world gait with digital technology? Validation, insights and recommendations from the Mobilise-D consortium.

[5] Palmerini, L., et al. "Mobility recorded by wearable devices and gold standards: the Mobilise-D procedure for data standardization." Scientific Data (2022).}

Notes

Todo-List

Prepare work:

• Linked this issue in the meta issue for algorithm development
• Created new branch and PR for this algorithm
• Input Output datatypes identified
• Relevant algorithm parameters identified
• Identified all relevant information
• Confirmed the info and general implementation plan with the rest of the developer team

Generate results from original implementation

• Original algorithm is running locally
• Example outputs from original algorithms generated and added to the repo

Implementation of Algorithm

• Core algorithm implemented
• Wrapper class implemented
• Initial peer-code review completed
• Extensive docstring added to algorithm class (should include references!)
• General tests added for the wrapper class
• Unit tests added (in particular to test edge cases)
• Example created (example should show how to use the algorithm, visualize the algo results and show a comparison the reference data and the results from the original implementation)
• Regression test for example created
• Added wrapper class to documentation
• Final code review completed
• PR merged and meta issue updated

This issue is meant to track the implementation of the {TD_ElGohary}

Short description: The turning detection toolbox contains code (MATLAB, R2018b) for detection of turning using angular accelerations recording with triaxial gyroscope worn/fixed on the lower back (close to body centre of mass)

project contains the code of the turning detection algorithm decribed in El-Gohary et al. 2014.

Four output files will be stored in the outdir in .json and .mat format, respectively, and both with an algorithm generic and algorithm specific file name. The Turning_ElGohary algorithm provides as output the regions of turnings with start and end in seconds since the start of the recording / trial, and several related parameters like the turning angle or the peak angular velocity.

Reference papers
El-Gohary, M., Pearson, S., McNames, J., Mancini, M., Horak, F., Mellone, S., & Chiari, L. (2014). Continuous monitoring of turning in patients with movement disability. Sensors (Basel, Switzerland), 14(1), 356–369. https://doi.org/10.3390/s140100356

Bascially, should we port what we implemented here:

https://gitlab.insight-centre.org/mobilise-d/dmo_cleaning

The example data from HA002 seems to have extremly high values for the initial contact in some cases. This was spotted for test11. Needs further investigation

At the moment, we don't have any way to handle transformers that change kwargs somehow.
For example the resampler would need to change the sampling_rate_hz in the kwargs for all subsequent steps.

This issue is meant to track the implementation of the GSD_A - GSDIluz

Short description: The algorithm convolves the raw IMU signal with a sin-wave template to find potential regions that contain regular motion of pre-defined frequency ranges. Afterwards, multiple post-processing steps are applied to check the plausibility of the identified GSDs
Link to original implementation: https://github.com/mobilise-d/Mobilise-D-TVS-Recommended-Algorithms/tree/master/GSDA
Reference Papers: https://jneuroengrehab.biomedcentral.com/articles/10.1186/1743-0003-11-48

Notes

• The algorithm shifts the signal after filtering. (see here) I don't get why this is needed. Is this to compensate for filter artifacts at the beginning of the signal?
• The algorithm uses a ORDER 200 filter! That seems high and will likely produce a bunch of signal artifacts.

Todo-List

Prepare work:

• Linked this issue in the meta issue for algorithm development
• Created new branch and PR for this algorithm (#22)
• Input Output datatypes identified
  • The algorithm just need the raw IMU data (specifically only acc) and the sampling rate.
    However, the algorithms seems to be designed for 100 Hz, but I don't see a reason, why it should not work with other sampling rates
• Relevant algorithm parameters identified
  • The algorithm has a couple of main parameters and a bunch of filter parameters. For the latter I am not sure how we can make them easily configurable, as it might be really complicated to explain the effect of changing them.
  • Main Parameters: window_size_s, window_overlap, activity_thres, max_number_steps_per_s, min_number_steps_per_s, min_step_peak_height
  • Initial bandpass filter paras for activity recognition: filter_order, cutoff_low, cutoff_high
  • Postprocessing was upright: upright_dc_threshold,
  • Postprocessing no sit to stand: check_window_duration, allowed_v_acc_difference
  • Postprocessing min-bout length: min_length
• Identified all relevant information
• Confirmed the info and general implementation plan with the rest of the developer team

Generate results from original implementation

• Original algorithm is running locally
• Example outputs from original algorithms generated and added to the repo

Implementation of Algorithm

• Core algorithm implemented
• Wrapper class implemented
• Initial peer-code review completed
• Extensive docstring added to algorithm class (should include references!)
• General tests added for the wrapper class
• Unit tests added (in particular to test edge cases)
• Example created (example should show how to use the algorithm, visualize the algo results and show a comparison the reference data and the results from the original implementation)
• Regression test for example created
• Added wrapper class to documentation
• Final code review completed
• PR merged and meta issue updated

During the process of reimplementing the WBA, I removed all rules that are not relevant anymore considering that we don't filter WBs based on their elevation change or the existance of turns anymore.

In this context we also removed the concepts of events from the WBA to make things simpler. This can be reintroduced once the below rules are reimplemented.

However, It might still be interesting to implement versions of these rules in the future.
Below is the source code of these rules.
Note, that this code will not work with the version of the WBA implemented in gaitlink at the moment.
But they can serve as a starting point to finalize the implementations of these rules:

class EventCriteria:
    event_type: str

    def filter_events(self, event_list):
        if event_list is None:
            raise ValueError("You are using an event based Criteria, without providing any events.")
        events = next(
            (event["events"] for event in event_list if event["name"] == self.event_type),
            [],
        )
        events = self._filter_events(events)
        return events

    @staticmethod
    def _convert_to_start_stop(events):
        event_start_end = np.array([[event["start"], event["end"]] for event in events])
        event_start_end = event_start_end[event_start_end[:, 0].argsort()]
        return event_start_end

    def _filter_events(self, event_list: list[dict]) -> list[dict]:
        return event_list


# TODO: "Precompile relevant event list for performance
class EventTerminationCriteria(BaseWBCriteria, EventCriteria):
    """Terminate/prevent starts of WBs in case strides overlap with events.

    Parameters
    ----------
    event_type
        Which type of event (from the event list) should be used
    termination_mode
        This controls under which circumstances a stride should be considered invalid for the WB.

        "start": A stride is invalid, if since the beginning of the last stride a new event was started.
        "end": A stride is invalid, if an event was ended since the beginning of the last stride
        "both": A stride is invalid, if any of the above conditions applied
        "ongoing": A stride is invalid, if it has any overlap with any event.

        At the very beginning of a recording, the start of the last stride is equal to the start of the recording.

    """

    termination_mode: str

    _termination_modes = (
        "start",
        "end",
        "both",
        "ongoing",
    )

    def __init__(
        self,
        event_type: str,
        termination_mode: Literal["start", "end", "both", "ongoing"] = "ongoing",
        comment: Optional[str] = None,
    ) -> None:
        self.event_type = event_type
        self.termination_mode = termination_mode
        super().__init__(comment=comment)

    def check_wb_start_end(
        self,
        stride_list: list[dict],
        original_start: int,
        current_start: int,
        current_end: int,
        event_list: Optional[list[dict]] = None,
    ) -> tuple[Optional[int], Optional[int]]:
        last_stride_start = 0
        if current_end >= 1:
            last_stride_start = stride_list[current_end - 1]["start"]
        current_stride_end = stride_list[current_end]["end"]
        current_stride_valid = self._check_stride(last_stride_start, current_stride_end, event_list)
        if current_stride_valid:
            return None, None
        # If the current stride is not valid, we either need to delay the start, if no proper WB has been started yet,
        # or we need to terminate the WB.
        if current_end == current_start:
            # Prevent the start
            return current_end + 1, None
        # Terminate the WB
        return None, current_end - 1

    def _check_stride(self, last_stride_start: float, current_stride_end: float, event_list):
        if self.termination_mode not in self._termination_modes:
            # We do this check here to avoid computing the events if the termination mode is invalid
            raise ValueError(f'"termination_mode" must be one of {self._termination_modes}')
        events = self.filter_events(event_list)
        if not events:
            return True
        event_start_end = self._convert_to_start_stop(events)

        events_started_since_last_stride = event_start_end[
            np.nonzero((event_start_end[:, 0] >= last_stride_start) & (event_start_end[:, 0] < current_stride_end))
        ]

        events_ended_since_last_stride = event_start_end[
            np.nonzero((event_start_end[:, 1] >= last_stride_start) & (event_start_end[:, 1] < current_stride_end))
        ]

        if self.termination_mode == "start":
            return len(events_started_since_last_stride) == 0
        if self.termination_mode == "end":
            return len(events_ended_since_last_stride) == 0
        if self.termination_mode == "both":
            return len(events_started_since_last_stride) == 0 and len(events_ended_since_last_stride) == 0
        if self.termination_mode == "ongoing":
            # Find events that where started before and are still ongoing
            ongoing_events = event_start_end[
                np.nonzero((event_start_end[:, 0] <= last_stride_start) & (event_start_end[:, 1] >= current_stride_end))
            ]
            return len(events_started_since_last_stride) == 0 and len(ongoing_events) == 0
        # We never reach this point
        raise ValueError()


class EventInclusionCriteria(BaseWBCriteria, EventCriteria):
    """Test if a WB is fully or partially covered by an event."""

    event_type: str
    overlap: str

    _overlap_types = ("partial", "contains", "is_contained", "no_overlap")

    def __init__(
        self,
        event_type: str,
        overlap: str = "partial",
        comment: Optional[str] = None,
    ) -> None:
        self.overlap = overlap
        self.event_type = event_type

        super().__init__(comment=comment)

    def check_include(self, wb: dict, event_list: Optional[list[dict]] = None) -> bool:
        # TODO: TEST
        events = self.filter_events(event_list)
        if not events:
            return True
        event_start_end = self._convert_to_start_stop(events)

        min_ends = np.minimum(event_start_end[:, 1], wb["end"])
        max_start = np.maximum(event_start_end[:, 0], wb["start"])
        amount_overlap = min_ends - max_start

        if self.overlap == "contains":
            return len(event_start_end[amount_overlap >= event_start_end[:, 1] - event_start_end[:, 0]]) > 0
        if self.overlap == "is_contained":
            return len(event_start_end[amount_overlap >= wb["end"] - wb["start"]]) > 0
        if self.overlap == "no_overlap":
            return len(event_start_end[amount_overlap > 0]) == 0
        if self.overlap == "partial":
            return len(event_start_end[amount_overlap > 0]) > 0
        raise ValueError(f'"overlap" must be one of {self._overlap_types}')
        
class LevelWalkingCriteria(BaseWBCriteria):
    """Test if WB has no more than N consecutive non-level strides.

    A WB is terminated if there are more than N consecutive strides that are not level walking.
    """
    max_non_level_strides: Optional[int]
    max_non_level_strides_left: Optional[int]
    max_non_level_strides_right: Optional[int]
    level_walking_threshold: float

    @property
    def _max_lag(self) -> int:
        non_none_vals = [
            val
            for val in [self.max_non_level_strides, self.max_non_level_strides_left, self.max_non_level_strides_right]
            if val
        ]
        return max(non_none_vals)

    def __init__(
        self,
        level_walking_threshold: float,
        max_non_level_strides: Optional[int] = None,
        max_non_level_strides_left: Optional[int] = None,
        max_non_level_strides_right: Optional[float] = None,
        field_name: str = "elevation",
        comment: Optional[str] = None,
    ) -> None:
        self.max_non_level_strides = max_non_level_strides
        self.max_non_level_strides_left = max_non_level_strides_left
        self.max_non_level_strides_right = max_non_level_strides_right
        if level_walking_threshold < 0:
            raise ValueError("`level_walking_threshold` must be >0.")
        self.level_walking_threshold = level_walking_threshold
        self.field_name = field_name
        super().__init__(comment=comment)

    def check_wb_start_end(
        self,
        stride_list: list[dict],
        original_start: int,
        current_start: int,
        current_end: int,
        event_list: Optional[list[dict]] = None,
    ) -> tuple[Optional[int], Optional[int]]:
        past_strides = stride_list[original_start : current_end + 1]
        consecutive_section = []
        for stride in reversed(past_strides):
            # We consider nan values always as level  walking!
            stride_height_change = abs(stride["parameter"][self.field_name])
            if not np.isnan(stride_height_change) and stride_height_change >= self.level_walking_threshold:
                consecutive_section.insert(0, stride)
            else:
                break
        if not consecutive_section:
            return None, None

        is_non_level = self._check_subsequence(consecutive_section)
        if is_non_level:
            # If we are at the beginning of the WB, we will change the start.
            if len(consecutive_section) == len(past_strides):
                return current_end + 1, None

            # If we are in the middle of a WB, we want to terminate it
            return None, current_end - len(consecutive_section)

        return None, None

    def _check_subsequence(self, stride_list) -> bool:
        """Check if the detected part exceeds our thresholds."""
        if self.max_non_level_strides is not None:
            return len(stride_list) >= self.max_non_level_strides
        if self.max_non_level_strides_left is None and self.max_non_level_strides_right is None:
            return False
        foot = [s["foot"] for s in stride_list]
        foot_count = Counter(foot)
        if self.max_non_level_strides_left is not None:
            return foot_count["left"] >= self.max_non_level_strides_left
        if self.max_non_level_strides_right is not None:
            return foot_count["right"] >= self.max_non_level_strides_right
        return False


class TurnAngleCriteria(EventTerminationCriteria):
    event_type: str = "turn"
    min_turn_angle: float
    max_turn_angle: float
    min_turn_rate: float
    max_turn_rate: float

    _serializable_paras = (
        "min_turn_angle",
        "max_turn_angle",
        "min_turn_rate",
        "max_turn_rate",
    )

    _rule_type: str = "turn_event"

    def __init__(
        self,
        min_turn_angle: Optional[float] = None,
        max_turn_angle: Optional[float] = None,
        min_turn_rate: Optional[float] = None,
        max_turn_rate: Optional[float] = None,
        comment: Optional[str] = None,
    ) -> None:
        self.min_turn_angle = min_turn_angle
        self.max_turn_angle = max_turn_angle
        self.min_turn_rate = min_turn_rate
        self.max_turn_rate = max_turn_rate
        super().__init__(
            event_type=self.event_type,
            termination_mode="ongoing",
            comment=comment,
        )

    def _filter_events(self, event_list: list[dict]) -> list[dict]:
        min_turn_angle, max_turn_angle = check_thresholds(
            self.min_turn_angle, self.max_turn_angle, allow_both_none=True
        )
        min_turn_rate, max_turn_rate = check_thresholds(self.min_turn_rate, self.max_turn_rate, allow_both_none=True)
        valid_events = []
        for e in event_list:
            if (min_turn_angle <= np.abs(e["parameter"]["angle"]) <= max_turn_angle) and (
                ((e["end"] - e["start"]) > 0)
                and min_turn_rate <= np.abs(e["parameter"]["angle"]) / (e["end"] - e["start"]) <= max_turn_rate
            ):
                valid_events.append(e)
        return valid_events

This issue is meant to track the implementation of the ICD_A

Short description: Algorithm implemented on a pre-processed vertical acceleration signal recorded on lower back. This signal is first detrended and then low-pass filtered (FIR, fc=3.2 Hz). The resulting signal is numerically integrated (cumtrapz) and differentiated using a Gaussian continuous wavelet transformation (CWT, scale 9, gauss2). The initial contact (IC) events are identified as the positive maximal peaks between successive zero-crossings.
Link to original implementation: https://github.com/mobilise-d/Mobilise-D-TVS-Recommended-Algorithms/tree/master/ICDA
Reference Papers:

• https://www.sciencedirect.com/science/article/pii/S0966636212000707?via%3Dihub
• https://ieeexplore.ieee.org/document/9176281
• https://jneuroengrehab.biomedcentral.com/articles/10.1186/s12984-023-01198-5
• https://www.nature.com/articles/s41597-023-01930-9

Notes

Todo-List

Prepare work:

• Linked this issue in the meta issue for algorithm development
• Created new branch and PR for this algorithm
• Input Output datatypes identified
• Relevant algorithm parameters identified
• Identified all relevant information
• Confirmed the info and general implementation plan with the rest of the developer team

Generate results from original implementation

• Original algorithm is running locally
• Example outputs from original algorithms generated and added to the repo

Implementation of Algorithm

• Core algorithm implemented
• Wrapper class implemented
• Initial peer-code review completed
• Extensive docstring added to algorithm class (should include references!)
• General tests added for the wrapper class
• Unit tests added (in particular to test edge cases)
• Example created (example should show how to use the algorithm, visualize the algo results and show a comparison the reference data and the results from the original implementation)
• Regression test for example created
• Added wrapper class to documentation
• Final code review completed
• PR merged and meta issue updated

For the wrist sensor (INDIP system), 'LeftWrist' or 'RightWrist' is stored as position in the data.mat file depending on the handedness of the participant.

It would be favourable to be able to load the dataset as follows by generically indicating 'Wrist':

mobilised_dataset = GenericMobilisedDataset(paths_list = paths, 
                                            raw_data_sensor = 'INDIP',
                                            sensor_positions = ('Wrist',),
                                            test_level_names = ("timemeasure", "recording")
                                            )

We should add methods that simplify working with the typical input data.

THis should include functions that check the correct dtype, columns, etc.

Similar (but likely simpler) to https://github.com/mad-lab-fau/gaitmap/blob/master/gaitmap/utils/_datatype_validation_helper.py

We developed Mobilsed internal tooling to plot the output of the WBA in an interactive way using Bokeh. This was a nice was to inspect per stride results and spot potential issues with WBA rules or algorithms in general.

We might want to port this tooling over as well.

In the current pipeline, all algorithms output values per second.
A second block then interpolates these values to strides.

This is somewhat of a confusing process. Originally, we did it this way to support algorithms with different outputs (steps/strides/per second). We should reevaluate and maybe find a better way of handling that.

For the main pipeline we will still need values per stride. Otherwise the entire WBassembley stuff will not work. However, I think on an algorithm/package level, we should support that people might want to use different output formats.

One related issue is if we modify the WBA to be felxible enough to also handle just per-sec inputs. Depending on the rules that are used this should be possible, but I am not sure how much sense this makes.

If we want to make per-sec output a first-class citizen, we should consider this...

The name "Gaitlink" was used as place holder name for now.

We should have a proper discussion on the name of the new pipeline project in (https://github.com/mobilise-d/pipeline_rewrite_planning/issues/1)

Once resolved, we need to rename the package and the repo

Our cadence deteciton methods basically use IC detectors and then a ourlier correction (implemented in the matlab cad2sec function) to calculate Cadence.

cad2sec is not that easy to reimplement in matlab.

Some components we need:

def hampel_filter_vectorized(input_series: np.ndarray, window_size: int, n_sigmas: float = 3.0) -> np.ndarray:
    k = 1.4826  # Scale factor for Gaussian distribution
    new_series = input_series.copy()
    
    # Create the median filtered series
    median_series = median_filter(input_series, size=window_size, mode='reflect')
    # Calculate the median absolute deviation with the corrected function
    scaled_mad = k * median_filter(median_abs_deviation(input_series, scale='normal'), size=window_size, mode='reflect')

    # Detect outliers
    outliers = np.abs(input_series - median_series) > n_sigmas * scaled_mad
    
    # Replace outliers
    new_series[outliers] = median_series[outliers]

    return new_series

In _parse_single_sensor_data, acc data is converted by multiplying with 9.81 (m/s^2). If the gaitlink MobiliseD dataset is used with external algorithms which also perform a conversion, this automatic conversion by gaitlink might not be wished.

https://github.com/mobilise-d/gaitlink/blob/eb7091661c1cacb47deeca02d98101eb0a40f814/gaitlink/data/_mobilised_matlab_loader.py#L348C3-L348C3

@AKuederle : Can an additional flag be used to trigger whether conversion is performed (default: conversion is performed)?

We should use a small dataset (1 or 2 participants) as concrete example data that we will store in the repo.
The data will also contain Ground Truth data that we can use to check the correctness

This data will be used for regression tests and examples.

We should also produce the results of the Matlab algorithms on this data and store these results in the repo as well as comparison.

In matlab GS have one sample more. It seems like the way they are cut in Matlab, the same GS has one more sample at the beginning. We need to investigate why this is the case, and load the GS reference values correctly, so that the matlab and Python version match.

We need to implement the WP6 aggregation methods into Python for the CVS. We could try to directly implement them in a way that they are compatible with the new pipeline

Might be nice, to have post-processing steps as composable blocks as well.
For now this seems overkill.

For the future, here is what i implemented for testing, but hten removed again:

from typing import Optional

import pandas as pd
from tpcp import Algorithm
from typing_extensions import Self

from gaitlink.consts import GRAV
from gaitlink.data_transform import BaseFilter, FirFilter


class BaseGsdPostProcessor(Algorithm):
    _action_methods = ("post_process",)

    # Other paramters
    data: pd.DataFrame
    sampling_rate_hz: float

    # results
    processed_gsd_list_: pd.DataFrame

    def post_process(self, data: pd.DataFrame, gsd_list: pd.DataFrame, *, sampling_rate_hz: float, **kwargs) -> Self:
        raise NotImplementedError


class RemoveDuration(BaseGsdPostProcessor):
    def __init__(self, min_duration_s: Optional[float] = 5, max_duration_s: Optional[float] = None) -> None:
        self.min_duration_s = min_duration_s
        self.max_duration_s = max_duration_s

    def post_process(self, data: pd.DataFrame, gsd_list: pd.DataFrame, *, sampling_rate_hz: float, **_) -> Self:
        duration = (gsd_list["end"] - gsd_list["start"]) / sampling_rate_hz
        selected = pd.Series(True, index=gsd_list.index)
        if self.min_duration_s is not None:
            selected &= duration >= self.min_duration_s
        if self.max_duration_s is not None:
            selected &= duration <= self.max_duration_s
        self.processed_gsd_list_ = gsd_list[selected]
        return self


class RemoveNonUpright(BaseGsdPostProcessor):
    def __init__(
        self,
        upright_threshold_ms2: float = 0.5 * GRAV,
        dc_filter: BaseFilter = FirFilter(order=5, cutoff_freq_hz=0.5, filter_type="lowpass"),
    ) -> None:
        self.upright_threshold_ms2 = upright_threshold_ms2
        self.dc_filter = dc_filter

    def post_process(self, data: pd.DataFrame, gsd_list: pd.DataFrame, *, sampling_rate_hz: float, **kwargs) -> Self:
        # We calculate the mean of the vertical acc component for each gsd
        # If the mean is above the threshold, we keep the gsd
        acc_v = self.dc_filter.clone().filter(data["acc_z"], sampling_rate_hz=sampling_rate_hz).transformed_data_
        mean_acc = [acc_v.iloc[start:end].mean() for start, end in gsd_list[["start", "end"]].itertuples(index=False)]
        selected = pd.Series(mean_acc, index=gsd_list.index) > self.upright_threshold_ms2

        self.processed_gsd_list_ = gsd_list[selected]
        return self


class RemoveTransitions(BaseGsdPostProcessor):
    def __init__(self, time_window_s: float = 1, allowed_difference_per: float = 15) -> None:
        self.time_window_s = time_window_s
        self.allowed_difference_per = allowed_difference_per

    def post_process(self, data: pd.DataFrame, gsd_list: pd.DataFrame, *, sampling_rate_hz: float, **kwargs) -> Self:
        # We calculate the mean of the vertical acc component for the first and the last n seconds of each gsd
        # If the difference is above the threshold, we remove the gsd
        acc_v = data["acc_z"]
        for start, end in gsd_list[["start", "end"]].itertuples(index=False):
            acc_v.iloc[start:end]
            
        # TODO: Finsih that

This issue is meant to track the implementation of the {LR_Det_Ullrich} algorithm for left/right foot detection for a single IMU sensor placed on the lower back.

Short description:

1) The McCamley algorithm

In the original McCamley algorithm, the angular velocity around the vertical axis ($gyr_{v}$) serves as the distinguishing factor for identifying left and right ICs. The process involves the following steps:

• Signal Pre-processing: Subtracting the signal mean and applying a low-pass filter (4th order Butterworth filter with a 2 Hz cut-off frequency).
• IC Assignment: Analyzing the sign of the filtered $gyr_{v}$ value at the IC time point $n$ for classification. If the value is positive, the IC is attributed to the right foot; if negative, it's attributed to the left foot.

As a first extension to the original McCamley algorithm, the angular velocity around the anterior-posterior axis, $gyr_{ap}$, can resemble a periodic wave with a constant phase shift w.r.t. $gyr_{v}$ after application of the low-pass filter described above. This is also a suitable input signal for the McCamley algorithm, when inverting the sign. A second and final extension to the original McCamley algorithm is to use the combination of the filtered signals for the vertical and anterior-posterior signals, $gyr_{comb}$:

$$ \begin{equation} gyr_{comb} = gyr_{v} - gyr_{ap} \end{equation} $$

2) Machine Learning Approaches

In the ML-based algorithms, we expand the feature set by incorporating the first and second derivatives of the filtered signals at the time points of the $n$ ICs. Consequently, for a dataset containing a total of $N$ ICs, this results in a $N \times 6$ feature matrix. To ensure uniformity, the feature set is min-max normalized.

Four different algorithms are employed for the left/right detection, which is essentially a binary classification task:

• Linear Support Vector Machine (SVM-lin)
• Radial Basis Function Support Vector Machine (SVM-rbf)
• k-Nearest Neighbours (kNN)
• Random Forest Classifiers (RFC)

Link to original implementation: N/A

Reference Papers: Ullrich M, Kuderle A, Reggi L, Cereatti A, Eskofier BM, Kluge F. Machine learning-based distinction of left and right foot contacts in lower back inertial sensor gait data. Annu Int Conf IEEE Eng Med Biol Soc. 2021, available at: https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=9630653

Notes

TBA

Todo-List

Provide ground truth labels for each IC (left or right) from the gold-standard data. Use this to compute the accuracy of the algorithms on the test data. The accuracy is formally defined as the number of agreements between predicted left and right labels and the ground truth labels divided by the total number of ICs, as described in below:

Prepare work:

• Linked this issue in the meta issue for algorithm development
• Created new branch and PR for this algorithm
• Input Output datatypes identified
• Relevant algorithm parameters identified
• Identified all relevant information
• Confirmed the info and general implementation plan with the rest of the developer team

Generate results from original implementation

• Original algorithm is running locally
• Example outputs from original algorithms generated and added to the repo

Implementation of Algorithm

• Core algorithm implemented
• Wrapper class implemented
• Initial peer-code review completed
• Extensive docstring added to algorithm class (should include references!)
• General tests added for the wrapper class
• Unit tests added (in particular to test edge cases)
• Example created (example should show how to use the algorithm, visualize the algo results and show a comparison the reference data and the results from the original implementation)
• Regression test for example created
• Added wrapper class to documentation
• Final code review completed
• PR merged and meta issue updated

In multiple cases we have predefined sets of configs for a specific algorithm. E.g. parameters specified on a specific dataset, or algorithm parameters proposed for a specific cohort.

In these cases, we need to provide these pre-defined configs to the user somehow.
There are a couple of options on how to do that. I will try to explain from the user perspective how these might look like.

Let's assume we have a GSD algorithm (GsdSomething), which takes a 2 parameters para_1 and para_2. They are both set to sensible defaults (para_1=1, para_2=10), but we want to ship different configurations of these parameters with gaitlink.
Specifically, we know that for HA the following configuration would work best para_1=3, para_2=7
and for MS patients the following config would work best: para_1=1, para_2=9.

So the question is how do we provide these options to the user.

1. As a Enum of config options. This enum of config options could be provided as a class var on the algorithm itself (or seperatly, but I think I prefer on the object)

from gaitlink.gsd import GsdSomething

# Initialize the default
gsd = GsdSomething()

# Initialize the HA optimized
gsd = GsdSomething(**GsdSomething.PRE_TRAINED.ha)

# Initialize the MS optimized
gsd = GsdSomething(**GsdSomething.PRE_TRAINED.ms)

# Take the HA default as a starting point and modify some parameters
# 1. Option
config = GsdSomething.PRE_TRAINED.ha
config["para_1"] = 3
gsd = GsdSomething(**config)

# 2. Option
gsd = GsdSomething(**GsdSomething.PRE_TRAINED.ms)
gsd = gsd.set_params(para_1 = 3)

2. As pre-initialized instances

Basically we store pre-initialized objects with the different parameters instead of the parameters.
One step less for the user, but could get complicated, when initializing the object has some runtime cost (e.g. large ML model)

from gaitlink.gsd import GsdSomething

# Initialize the default
gsd = GsdSomething()

# Get HA optimized (this is directly the model object)
gsd = GsdSomething.PRE_TRAINED.ha

# Get MS optimized (this is directly the model object)
gsd = GsdSomething.PRE_TRAINED.ms

# Take the HA default as a starting point and modify some parameters
# Only the second option still works
gsd = GsdSomething.PRE_TRAINED.ha
gsd = gsd.set_params(para_1 = 3)

3. Using some loader function

Basically the same as 2., but instead of pre-initializing we have a class-constructure method that we can pass the info about which version we want to load.

from gaitlink.gsd import GsdSomething

# Initialize the default
gsd = GsdSomething()

# Get HA optimized (this is directly the model object)
gsd = GsdSomething.from_pretrained("ha")

# Get MS optimized (this is directly the model object)
gsd = GsdSomething.from_pretrained("ms")

# Take the HA default as a starting point and modify some parameters
# Only the second option still works
gsd = GsdSomething.from_pretrained("ha")
gsd = gsd.set_params(para_1 = 3)

I am not sure what approach would be the nicest from a user perspective. Option 1 and 2 have the advantage, that we could implement them in a way, where you can get autocomplete of the model names in the IDE.
Option 2 seems to be the least explizit, and hence, likely most confusing.
The big advantage of option 3 would be that we could easily load configs or models from files without doing some property magic in the background. But we should decide from a user-perspective IMO.

I am personally leaning towards option 1, but not 100 % sure.
Any opinions?

Also @felixkluge @rmndrs89

Following the Left/Right Wrist implementation, I encountered the following issue when trying to load the complete TVS dataset. Currently, loading seems to be too strict when specifying missing_sensor_error_type=ignore/warn. E.g., if one data point does not contain a sensor (e.g., specified sensor is not available due to technical issues during data collection), an error is raised, hindering processing of the rest of the dataset. Expected behaviour would rather be to ignore this error, setting imu_data = {}, and to continue with the next data point.

Example case:

• Creating Dataset with sensor_positions=("LeftWrist", "RightWrist"),
• missing_sensor_error_type=ignore/warn
• A data point in the data set contains neither LeftWrist nor RightWrist
  -> Then, an error is raised, because all_imu_data == {}

This concerns the following lines:
https://github.com/mobilise-d/gaitlink/blob/362b14d72e04c77b32f5f65665fbd0f282c7c53a/gaitlink/data/_mobilised_matlab_loader.py#L289C9-L293C1

Suggestions:
In case missing_sensor_error_type=ignore/warn is set, no error should be raised, even when no imu data is available

We could replace the above line by:

if all_imu_data == {}:
            error_message = (
                f"Expected at least one valid sensor position for {raw_data_sensor}. Given: {sensor_positions}"
            )
            if missing_sensor_error_type == "raise":
                raise ValueError(error_message)
            if missing_sensor_error_type == "warn":
                warnings.warn(error_message, stacklevel=1)

The only issue is that, in case no other sensor (position) has been specified, that the set of sampling_rate_valueswill be empty, leading to the Expected all sensors across all positions to have the same sampling rate, but found {sampling_rates}error. Here, a check whether any sampling rate has been specified is also needed.

BaseGsdDetector -> BaseGsDetector
gsd_list -> gs_list

The d stands for detector

Follow-up on #69 to fix the review comments

At the moment we have gaitmap as a dependcy to easily pull utility functions over.

However, we should remove gaitmap dependencies at some point and port the relevant functions over to this repo

This issue should serve as a place for general Q&A until we find a better place to document these topics.

Feel free to ask new questions below. Curated answers will be provided by editing the top post.

Q&A

Why per-second value output for CAD/SL/Gaitspeed?

While we only implement Physics based CAD and SL algorithms, which can provide step-level output, in the original algorithm comparison, we also used ML algorithm that estimated the respective values per window of data. The only common ground was to interpolate step based values to per-second values.
This seems a little awkward, and also introduces some level of error, but allows to easily replace the physics based algorithms with ML based algorithms.

WB vs LWB vs GS

Fundamentally, we differentiate between walking bouts and gait sequences. Gait sequences are "regions likely to contain gait". They don't follow strict rules, it is basically "whatever the GS detection algorithm think its gait". Walking bouts follow strict defined rules. These rules were defined as part of a consensus process within Mobilise-D and cover things like the minimum number of strides, the allowed break between two strides and others. Because these rules can only be applied once all parameters are calculated, the final WB can be significantly different from the original GSs.
We further differentiate between Walking Bout and Level Walking Bout (LWB). The latter is a WB that contains not incline or decline walking. The idea behind this is, that clinically comparable gait parameters can likely only extracted from level walking. However, because the detection of inclines and declines from just a single lower back IMU is challenging, we rarely use this definition and usually do all comparisons on the WB.

Historic Naming: In older standardized data, WBs were referred to as "Continuous Walking Period (CWP)" or "MacroWbs" and Level WBs were called "MicroWBs".

Why do all algorithms store the input on the object? Does this increase memory consumption?

All algorithm objects make the input data available via self.data and self.ic_list (and others) available after calling the action method.
This is a "convention" that leads to some nice side-effects. This way, the final object has all the information about inputs and outputs. You could write for example a plot func, that just takes the algo instance as input and have all the information available that you need. I also thought about memory consumption in this context, and the reality is, that it won't matter in 99% of the cases. Just storing the object on the instance will not result in an increase in memory consumption. It only stores the reference to the object. This means the only thing that "could" cause issues here, is that as long as the instance with the results attached exists, the original data will not be cleared from memory. But, it basically never happens, that the original data is out of scope and an algorithm object is still in scope, blocking the data deletion. I played around with that a bit (as we are doing the same in gaitmap) and it never was an issue in any of our usecases.

The Picture above shows the full pipeline we need to implement. For most algorithmic steps this means one algorithm. Only for GSD and CAD there are two.

In addition to the blocks shown above, there is a Reorientation Block that should catch instances where the signal is wrongly oriented.

Overview

• Reorientation block (#112 )
• GSD 1 (#21)
• GSD 2 (#24)
• Turning (#31)
• IC-Detection (#28 )
• CAD 1/IC 2 (#29 #72 )
• CAD 2/IC 3 (#29 #84 )
• Cad2Sec (#61)
• #115
• Left Right - Normal (#32, #100 )
• Left Right - ML pretrained (#32)
• Walking Speed ()
• Stride interpolation ()
• WBA (#35 )
• Final Parameter Formatting ()
• Dmo cleaning (#12, #91 )
• Aggregation Methods (#13)

How to implement new algorithms

When implementing a new algorithm, create a new issue with the template and update this issue with a link to all relevant issues and PRs.

1. Find the old implementation and try to understand it. What is the actual core algorithm? What is just pre-processing? What are the relevant parameters of the algorithm? What is the output type (per-second, per-step)? What inputs are required (raw data, sensor height, detected steps, ...)?
2. Get the old version to run offline on the example data. "Mock" previous outputs using the ground truth. Store the results somewhere
3. Discuss in the group, before starting any implementation. Are the outputs as expected? Do we understand how the algorithm is supposed to work?
4. Start building an Example (i.e. a notebook or Python file to prototype). Have your raw inputs and your example output from before there to play around.
5. Create a first implementation of the core functionality. Just build a simple Python function. Ignore class interface for now (I would suggest to try AI-tools to make the initial conversion from matlab to Python). Check if your implementation gets the same output as before. If there are deviations, that can not/should not be "fixed" document them
6. Build the class interface + tests + docu + example. Example should cover basic usage of algorithm + comparison to old matlab output. Write down anything that could interesting for someone in the future.

To avoid running subsequent (as in after the IC detection) on regions of signal that are not gait, we use the concept of "refined gait sequences" in Mobilise-D, where we cut a gait sequence to the region where valid ICs were detected. In the orignal pipeline that was an inherit output of the IC detectors.

We could at a repective property to the IC detection Baseclass (e.g. refined_gait_region_/refined_gait_region_start_end_), that would provide the new start ends and the cut data.

This could then be used in subsequent processing steps.
The difficult part will be to manage index offset in the subsequent processing steps.

Absolutely not necessary, but fun.

Creative ideas welcome :)

This issue is meant to track the implementation of the {GSD_LowBackAcc}

Short description: {The Gait Sequence Detection toolbox contains code (MATLAB, R2018b) for detection of gait (walking) sequences using body acceleration recorded with a triaxial accelerometer worn/fixed on the lower back (close to body center of mass).

The algorithm was developed and validated using data recorded in patients with impaired mobility (Parkinson’s disease, multiple sclerosis, hip fracture, post-stroke and cerebral palsy).

The algorithm detects the gait sequences based on identified steps. First, the norm of triaxial acceleration signal is detrended and low-pass filtered (FIR, fc=3.2Hz). In order to enhance the step-related features (peaks in acceleration signal) the obtained signal is further processed using continuous wavelet transform, Savitzky-Golay filters and Gaussian-weighted moving average filters [2]. The ‘active’ periods, potentially corresponding to locomotion, are roughly detected and the statistical distribution of the amplitude of the peaks in these active periods is used to derive an adaptive (data-driven) threshold for detection of step-related peaks. Consecutive steps are associated to gait sequences [1, 2].

The gait sequence detection algorithm is implemented in the main function GSD_LowerBackAcc.m. The script example_GSD.m contains an exemplary application of the algorithm, using input data in the standardized format adopted by MobiliseD project [4]. However, the algorithm can be applied for any data in the specified format (see Input/Output description).

Note that this algorithm is referred as GSDB in the validation study [3].}
Link to original implementation: {Link}
Reference Papers: {[1] Paraschiv-Ionescu, A, et al. "Locomotion and cadence detection using a single trunk-fixed accelerometer: validity for children with cerebral palsy in daily life-like conditions." Journal of neuroengineering and rehabilitation 16.1 (2019): 1-11.

[2] Paraschiv-Ionescu, A, Soltani A, and Aminian K. "Real-world speed estimation using single trunk IMU: methodological challenges for impaired gait patterns." 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). IEEE, 2020.

[3] Micó-Amigo, M. E., Bonci, T., Paraschiv-Ionescu, A., Ullrich, M., Kirk, C., Soltani, A., ... & Del Din, S. (2022). Assessing real-world gait with digital technology? Validation, insights and recommendations from the Mobilise-D consortium.

[4] Palmerini, L., et al. "Mobility recorded by wearable devices and gold standards: the Mobilise-D procedure for data standardization." Scientific Data (2022).}

Notes

Todo-List

Prepare work:

• Linked this issue in the meta issue for algorithm development
• Created new branch and PR for this algorithm
• Input Output datatypes identified
• Relevant algorithm parameters identified
• Identified all relevant information
• Confirmed the info and general implementation plan with the rest of the developer team

Generate results from original implementation

• Original algorithm is running locally
• Example outputs from original algorithms generated and added to the repo

Implementation of Algorithm

• Core algorithm implemented
• Wrapper class implemented
• Initial peer-code review completed
• Extensive docstring added to algorithm class (should include references!)
• General tests added for the wrapper class
• Unit tests added (in particular to test edge cases)
• Example created (example should show how to use the algorithm, visualize the algo results and show a comparison the reference data and the results from the original implementation)
• Regression test for example created
• Added wrapper class to documentation
• Final code review completed
• PR merged and meta issue updated