{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":41880,"databundleVersionId":5677426,"sourceType":"competition"}],"dockerImageVersionId":30627,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 0. Preamble","metadata":{}},{"cell_type":"markdown","source":"## 0.1. Imports","metadata":{}},{"cell_type":"code","source":"# utils\nimport os, gc\nimport numpy as np\nimport pandas as pd \nfrom functools import partial, update_wrapper\nfrom tqdm.auto import tqdm\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\nfrom graphviz import Digraph\nimport gc\n\n# sklearn\nfrom sklearn.model_selection import train_test_split, StratifiedGroupKFold\nfrom sklearn.metrics import accuracy_score, average_precision_score\nfrom sklearn.preprocessing import StandardScaler, RobustScaler\nfrom sklearn.decomposition import FastICA\n\n# tensorflow\nfrom tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau, CSVLogger\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom keras import layers\nfrom keras.preprocessing.sequence import pad_sequences\n\n# scipy\nfrom scipy.signal import welch, spectrogram, butter, filtfilt, hilbert\nfrom scipy import signal","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-02-25T11:48:40.144367Z","iopub.execute_input":"2024-02-25T11:48:40.144779Z","iopub.status.idle":"2024-02-25T11:48:58.576963Z","shell.execute_reply.started":"2024-02-25T11:48:40.144743Z","shell.execute_reply":"2024-02-25T11:48:58.575620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 0.2. Configs","metadata":{}},{"cell_type":"code","source":"class CFG:\n    # debug mode\n    debug = False\n    \n    # parallel\n    num_workers = -1\n    \n    # dataset\n    competition_name = \"tlvmc-parkinsons-freezing-gait-prediction\"\n    target_cols = [\"StartHesitation\", \"Turn\", \"Walking\"]\n    seq_len = 5000\n    shift = 2500\n    offset = 1250\n    input_dir = Path(\"../input/\") / competition_name\n    output_dir = Path(\"./\")","metadata":{"execution":{"iopub.status.busy":"2024-02-25T11:48:58.579060Z","iopub.execute_input":"2024-02-25T11:48:58.580103Z","iopub.status.idle":"2024-02-25T11:48:58.588435Z","shell.execute_reply.started":"2024-02-25T11:48:58.580067Z","shell.execute_reply":"2024-02-25T11:48:58.586912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 0.3. Utils","metadata":{}},{"cell_type":"code","source":"def get_filename_no_ext(file_path):\n    if isinstance(file_path, str):\n        base = os.path.basename(file_path)\n    elif isinstance(file_path, Path):\n        base = file_path.name\n    else:\n        raise ValueError(\"Unsupported file path type. Please provide a string or PosixPath object.\")\n    \n    filename = os.path.splitext(base)[0]\n    return filename\n\n\ndef _cartesian_to_spherical(x, y, z):\n    '''\n    Given a 3D vector in cartesian coordinates, convert to spherical coordinates\n\n    Args:\n        x, y, z: float\n            x, y, z coordinates of the vector\n    \n    Returns:\n        r, theta, phi: float\n    '''\n    r = np.sqrt(x**2 + y**2 + z**2)  # Radial distance\n    if r == 0:\n        theta = 0  # Avoid division by zero\n    else:\n        theta = np.arccos(z / r)  # Polar angle\n    phi = np.arctan2(y, x)  # Azimuth angle\n    return r, theta, phi\n\n\ndef cartesian_to_spherical(data):\n    '''\n    Given a 2D array of (3, n_timepoints) in cartesian coordinates, \n    convert cartesian coordinates to spherical coordinates.\n\n    Args:\n        data: numpy array\n            2D array of (3, n_timepoints) in cartesian coordinates\n\n    Returns:\n        r, theta, phi: numpy array\n            2D array of (3, n_timepoints) in spherical coordinates\n    '''\n    r, theta, phi = np.apply_along_axis(lambda m: _cartesian_to_spherical(*m), \n                                        axis=0, arr=data)\n    return np.vstack((r, theta, phi))\n\n\ndef butter_bandpass_filter(data, lowcut, highcut, fs, order=5):\n    '''\n    Perform bandpass filtering on a 2D array of (n_channels, n_timepoints).\n\n    Args:\n        data: numpy array\n            2D array of (n_channels, n_timepoints)\n        lowcut: float\n            Lower cutoff frequency\n        highcut: float\n            Upper cutoff frequency\n        fs: float\n            Sampling frequency\n        order: int\n            Order of the filter\n    \n    Returns:\n        data: numpy array of (n_channels, n_timepoints)\n    '''\n    nyq = 0.5 * fs\n    low = lowcut / nyq\n    high = highcut / nyq\n    # Get the filter coefficients\n    b, a = signal.butter(order, [low, high], btype='band')\n    data = np.apply_along_axis(lambda m: signal.lfilter(b, a, m), axis=1, arr=data)\n    return data\n\n\ndef moving_average(data, window_size):\n    '''\n    Perform moving average on a 2D array of (n_channels, n_timepoints).\n\n    Args:\n        data: numpy array\n            2D array of (n_channels, n_timepoints)\n        window_size: int\n            Size of the moving average window\n\n    Returns:\n        data: numpy array of (n_channels, n_timepoints)\n    '''\n    kernel = np.ones(window_size) / window_size\n    data = np.apply_along_axis(lambda m: np.convolve(m, kernel, mode='same'), \n                               axis=1, arr=data)\n    return data\n\n\ndef window_sinc_convolve(data, window_size, fs):\n    '''\n    Perform windowed sinc convolution on a 2D array of (n_channels, n_timepoints).\n\n    Args:\n        data: numpy array\n            2D array of (n_channels, n_timepoints)\n        window_size: int\n            Size of the window\n        fs: float\n            Sampling frequency\n\n    Returns:\n        data: numpy array of (n_channels, n_timepoints)\n    '''\n    kernel = np.sinc(2 * window_size * np.arange(fs) / fs)\n    data = np.apply_along_axis(lambda m: np.convolve(m, kernel, mode='same'), \n                               axis=1, arr=data)\n    return data\n\n\ndef acceleration_to_velocity(data, fs):\n    '''\n    Given a 2D array of (n_channels, n_timepoints) in acceleration, \n    convert to velocity.\n\n    Args:\n        data: numpy array\n            2D array of (n_channels, n_timepoints) in acceleration\n        fs: float\n            Sampling frequency\n\n    Returns:\n        data: numpy array of (n_channels, n_timepoints) in velocity\n    '''\n    data = np.apply_along_axis(lambda m: np.cumsum(m) / fs, axis=1, arr=data)\n    return data\n\n\ndef normalize(data):\n    '''\n    Given a 2D array of (n_channels, n_timepoints), normalize each channel.\n\n    Args:\n        data: numpy array\n            2D array of (n_channels, n_timepoints)\n\n    Returns:\n        data: numpy array of (n_channels, n_timepoints)\n    '''\n    data = np.apply_along_axis(lambda m: (m - np.mean(m)) / np.std(m), axis=1, arr=data)\n    return data\n\n\ndef hilbert_transform(data):\n    '''\n    Given a 2D array of (n_channels, n_timepoints), perform hilbert transform \n    on each channel and return 2 arrays of (n_channels, n_timepoints) for\n    amplitude and phase.\n\n    Args:\n        data: numpy array\n            2D array of (n_channels, n_timepoints)\n\n    Returns:\n        amplitude: numpy array of (n_channels, n_timepoints)\n        phase: numpy array of (n_channels, n_timepoints)\n    '''\n    hilbert_transformed = np.apply_along_axis(lambda m: hilbert(m), axis=1, arr=data)\n    amplitude = np.abs(hilbert_transformed)\n    phase = np.angle(hilbert_transformed)\n    return amplitude, phase\n\n\ndef spectrogram(data, fs, nperseg=128, noverlap=64, max_freq=50, min_freq=0):\n    '''\n    Given a 2D array of (n_channels, n_timepoints), perform STFT on each channel\n    and return an array of frequencies, an array of time points, and a 3D array \n    of (n_channels, n_freqs, n_timepoints).\n\n    Args:\n        data: numpy array\n            2D array of (n_channels, n_timepoints)\n        fs: float\n            Sampling frequency\n        nperseg: int\n            Length of each segment\n        noverlap: int\n            Number of points to overlap between segments\n        max_freq: float\n            Maximum frequency to return\n        min_freq: float\n            Minimum frequency to return\n\n    Returns:\n        frequencies: numpy array of (n_freqs,)\n        timepoints: numpy array of (n_timepoints,)\n        spectrogram: numpy array of (n_channels, n_freqs, n_timepoints)\n    '''\n    frequencies, timepoints, spectrogram = signal.spectrogram(data, fs=fs, \n                                                              nperseg=nperseg, \n                                                              noverlap=noverlap)\n    # Remove frequencies above max_freq and below min_freq\n    frequencies = frequencies[(frequencies <= max_freq) & (frequencies >= min_freq)]\n    spectrogram = spectrogram[:, (frequencies <= max_freq) & (frequencies >= min_freq), :]\n    return frequencies, timepoints, spectrogram\n\n\ndef noise_remove_ica(data, n_components=3):\n    '''\n    Given a 2D array of (n_channels, n_timepoints), perform ICA to remove noise.\n\n    Args:\n        data: numpy array\n            2D array of (n_channels, n_timepoints)\n        n_components: int\n            Number of components to keep\n\n    Returns:\n        data: numpy array of (n_channels, n_timepoints)\n    '''\n    ica = FastICA(n_components=n_components)\n    data = ica.fit_transform(data.T).T\n    return data\n\n\ndef generate_shorter_sequences(long_x, long_y, window_size, stride):\n    '''\n    Generates shorter sequences from longer input sequences.\n\n    Args:\n        long_x: A numpy array representing the long input sequences with shape (n_samples, n_features).\n        long_y: A numpy array representing the corresponding long target sequences with shape (n_samples, n_targets).\n        window_size: An integer specifying the size of the shorter sequences to be generated.\n        stride: An integer specifying the stride or step size along the long sequences.\n\n    Yields:\n        A generator that produces tuples of shorter input and target sequences.\n\n    Example:\n        # Generate shorter sequences with window size 5 and stride 2\n        long_input = np.array([[1, 2, 3, 4, 5],\n                               [6, 7, 8, 9, 10],\n                               [11, 12, 13, 14, 15]])\n        long_target = np.array([[0, 1],\n                                [1, 0],\n                                [0, 1]])\n        for short_x, short_y in generate_shorter_sequences(long_input, long_target, window_size=3, stride=2):\n            print(\"Short Input:\", short_x)\n            print(\"Short Target:\", short_y)\n    '''\n    seq_len = long_x.shape[0]\n    for i in range(0, long_x.shape[0] - window_size + 1, stride):\n        yield long_x[i:i+window_size, :], long_y[i:i+window_size, :]\n    if (seq_len % window_size) % stride > 0 :\n        yield long_x[-window_size:, :], long_y[-window_size:, :]\n    \n\ndef plot_spectrogram(X, fs, maxfreq=5, minfreq=0, nperseg=512, ax=None, figsize=(12, 4)):\n    \"\"\"\n    Plots the spectrogram of a signal.\n\n    Args:\n    - X (array_like): The input signal.\n    - fs (float): The sampling frequency of the signal.\n    - maxfreq (float, optional): Maximum frequency to display in the spectrogram.\n    Default is 5 Hz.\n    - minfreq (float, optional): Minimum frequency to display in the spectrogram.\n    Default is 0 Hz.\n    - nperseg (int, optional): Length of each segment used to compute the FFT.\n    Default is 512.\n    - ax (matplotlib.axes.Axes, optional): The Axes object to plot on. If None,\n    a new figure will be created.\n    - figsize (tuple, optional): Figure size (width, height) in inches.\n    Default is (12, 4).\n\n    Return: None\n    \"\"\"\n    if ax is None:\n        fig, ax = plt.subplots(figsize=figsize)\n    f, t, Sxx = spectrogram(X, fs)\n    ax.pcolormesh(t, f, Sxx, cmap='magma')\n    ax.set_ylabel('Frequency [Hz]')\n    ax.set_xlabel('Time [sec]')\n    ax.set_ylim(minfreq, maxfreq)\n    \n    \nclass WrappedPartial:\n    \"\"\"\n    A class representing a wrapped partial function.\n\n    Parameters:\n    - original_func (callable): The original function to create a partial function from.\n    - *args: Positional arguments to fix in the partial function.\n    - **kwargs: Keyword arguments to fix in the partial function.\n    \"\"\"\n\n    def __init__(self, original_func, *args, **kwargs):\n        \"\"\"\n        Initialize a WrappedPartial instance.\n\n        Parameters:\n        - original_func (callable): The original function to create a partial function from.\n        - *args: Positional arguments to fix in the partial function.\n        - **kwargs: Keyword arguments to fix in the partial function.\n        \"\"\"\n        partial_func = partial(original_func, *args, **kwargs)\n\n        # Preserve __doc__ and __name__ attributes\n        update_wrapper(partial_func, original_func)\n        self.partial_func = partial_func\n\n    def getfunc(self, new_name=None):\n        \"\"\"\n        Get the partial function.\n\n        Parameters:\n        - new_name (str, optional): If provided, set the __name__ attribute of the partial function to this value.\n\n        Returns:\n        - callable: The partial function.\n        \"\"\"\n        partial_func = self.partial_func\n        if new_name is not None:\n            partial_func.__name__ = new_name\n        return partial_func \n    \n    \n    \nclass FunctionNode:\n    \"\"\"\n    A class representing a node associated with a function.\n\n    Parameters:\n    - func (callable): The function associated with the root node.\n    - children (list, optional): List of child nodes. Defaults to an empty list.\n    \"\"\"\n\n    def __init__(self, func, children=None):\n        \"\"\"\n        Initialize a FunctionTree instance.\n\n        Parameters:\n        - func (callable): The function associated with the root node.\n        - children (list, optional): List of child nodes. Defaults to an empty list.\n        \"\"\"\n        self.func = func\n        self.children = children or []\n\n    def add_child(self, child_node):\n        \"\"\"\n        Add a child node to the root node.\n\n        Parameters:\n        - child_node (FunctionTree): The child node to be added.\n        \"\"\"\n        self.children.append(FunctionNode(child_node))\n\n    def find_node(self, target_func):\n        \"\"\"\n        Recursively search for a node with a specific function in the tree.\n\n        Parameters:\n        - target_func (callable): The target function to search for.\n\n        Returns:\n        - FunctionTree or None: The node with the target function, or None if not found.\n        \"\"\"\n        if self.func == target_func:\n            return self\n        else:\n            for child in self.children:\n                found_node = child.find_node(target_func)\n                if found_node:\n                    return found_node\n        return None\n\n    def add_child_to_node(self, target_func, new_child):\n        \"\"\"\n        Add a child node to an arbitrary node in the tree.\n\n        Parameters:\n        - target_func (callable): The function associated with the target node.\n        - new_child (FunctionTree): The child node to be added to the target node.\n        \"\"\"\n        target_node = self.find_node(target_func)\n        if target_node:\n            target_node.add_child(new_child)\n        else:\n            print(f\"Node with function {target_func} not found.\")\n\n    def _evaluate(self, input_value):\n        \"\"\"\n        Recursively evaluate the tree and return a list of outputs from all leaf nodes.\n\n        Parameters:\n        - input_value: The input value to be used in the function evaluations.\n\n        Returns:\n        - list: A list of outputs from all leaf nodes.\n        \"\"\"\n        if not self.children:\n            return self.func(input_value), \n        else:\n            curr_result = self.func(input_value)\n            child_results = [child.evaluate(curr_result) for child in self.children]\n            \n            return  [result for sublist in child_results for result in sublist]\n        \n    def evaluate(self, input_value):\n        \"\"\"\n        Evaluate the tree and return a list of outputs from all leaf nodes.\n\n        Parameters:\n        - input_value: The input value to be used in the function evaluations.\n\n        Returns:\n        - numpy.ndarray: A numpy array stacking all data produced from the leaf nodes on axis 0.\n        The resulting shape is (n_leaves, ...).\n        \"\"\"\n        results = self._evaluate(input_value)\n        return np.array(results)\n        \n\n    def visualize(self, graph=None, parent_name=None, graphviz=None, size=None):\n        \"\"\"\n        Generate a graphical representation of the tree using Graphviz.\n\n        Parameters:\n        - graph (Digraph, optional): The Graphviz graph. Defaults to None.\n        - parent_name (str, optional): The name of the parent node. Defaults to None.\n        - graphviz (Digraph, optional): The original Graphviz graph. Defaults to None.\n        - size (tuple, optional): The size of the output graph. Defaults to None.\n\n        Returns:\n        - Digraph: The Graphviz graph.\n        \"\"\"\n        if graph is None:\n            graph = Digraph(format='png')\n\n            # Set the size if provided\n            if size:\n                graph.attr(size=size)\n\n            graphviz = graph\n        current_name = str(id(self))\n        graph.node(current_name, label=str(self.func.__name__))\n\n        if parent_name is not None:\n            graph.edge(parent_name, current_name)\n\n        for i, child in enumerate(self.children):\n            child.visualize(graph, current_name, graphviz=graphviz, size=size)\n\n        return graph\n    \n    \n\n","metadata":{"execution":{"iopub.status.busy":"2024-02-25T11:49:02.080787Z","iopub.execute_input":"2024-02-25T11:49:02.081202Z","iopub.status.idle":"2024-02-25T11:49:02.137915Z","shell.execute_reply.started":"2024-02-25T11:49:02.081173Z","shell.execute_reply":"2024-02-25T11:49:02.136131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 0.4. Data","metadata":{}},{"cell_type":"code","source":"# Dataframes\n\nsample_submission = pd.read_csv(CFG.input_dir / \"sample_submission.csv\")\ntdcsfog_metadata = pd.read_csv(CFG.input_dir / \"tdcsfog_metadata.csv\")\ndefog_metadata = pd.read_csv(CFG.input_dir / \"defog_metadata.csv\")\ndaily_metadata = pd.read_csv(CFG.input_dir / \"daily_metadata.csv\")\nsubjects = pd.read_csv(CFG.input_dir / \"subjects.csv\")\nevents = pd.read_csv(CFG.input_dir / \"events.csv\")\ntasks = pd.read_csv(CFG.input_dir / \"tasks.csv\")\n\n# Directories\n\ntdcsfog_train_dir = CFG.input_dir / 'train/tdcsfog'\ndefog_train_dir = CFG.input_dir / 'train/defog'\nnotype_train_dir = CFG.input_dir / 'train/notype'\n\ntdcsfog_test_dir = CFG.input_dir / 'test/tdcsfog'\ndefog_test_dir = CFG.input_dir / 'test/defog'\n\n\n# Subjects for each dataset\ntdcsfog_train_sessions = [get_filename_no_ext(p) for p in tdcsfog_train_dir.glob('*.csv')]\ndefog_train_sessions = [get_filename_no_ext(p) for p in defog_train_dir.glob('*.csv')]\nnotype_train_sessions = [get_filename_no_ext(p) for p in notype_train_dir.glob('*.csv')]\n\ntdcsfog_test_sessions = [get_filename_no_ext(p) for p in tdcsfog_test_dir.glob('*.csv')]\ndefog_test_sessions = [get_filename_no_ext(p) for p in defog_test_dir.glob('*.csv')]","metadata":{"execution":{"iopub.status.busy":"2024-02-25T11:49:12.182850Z","iopub.execute_input":"2024-02-25T11:49:12.183285Z","iopub.status.idle":"2024-02-25T11:49:13.057590Z","shell.execute_reply.started":"2024-02-25T11:49:12.183252Z","shell.execute_reply":"2024-02-25T11:49:13.056707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission","metadata":{"execution":{"iopub.status.busy":"2024-02-25T11:49:22.737782Z","iopub.execute_input":"2024-02-25T11:49:22.738167Z","iopub.status.idle":"2024-02-25T11:49:22.767555Z","shell.execute_reply.started":"2024-02-25T11:49:22.738140Z","shell.execute_reply":"2024-02-25T11:49:22.766313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. Designing the reprocessing and feature engineering pipeline\n\n","metadata":{}},{"cell_type":"code","source":"def fog_pipeline(fs):\n    # define partial functions so that they only accepts arrays\n    bandpass_delta = WrappedPartial(butter_bandpass_filter, lowcut=0.5, highcut=4, fs=fs).getfunc('bandpass_delta')\n    bandpass_theta = WrappedPartial(butter_bandpass_filter, lowcut=4, highcut=8, fs=fs).getfunc('bandpass_theta')\n    bandpass_alpha = WrappedPartial(butter_bandpass_filter, lowcut=8, highcut=16, fs=fs).getfunc('bandpass_alpha')\n    average_20_hz = WrappedPartial(moving_average, window_size=fs//5).getfunc('average_20_hz')\n    acceleration = WrappedPartial(lambda x : x).getfunc('acceleration')\n    velocity = WrappedPartial(acceleration_to_velocity, fs=fs).getfunc('velocity')\n\n    # define tuple-handing leaf-node functions for hilbert\n    amplitude = lambda htup: htup[0] \n    phase = lambda htup: htup[1] \n\n    # Constructing a function tree, starting with normalization\n    preproc_tree = FunctionNode(normalize)\n\n    # Level 1 functions\n    preproc_tree.add_child(bandpass_delta)\n    preproc_tree.add_child(bandpass_theta)\n    preproc_tree.add_child(bandpass_alpha)\n    preproc_tree.add_child(average_20_hz)\n\n\n    # Level 2 functions\n    for band, child in zip(['delta', 'theta', 'alpha'],\n                           [bandpass_delta, bandpass_theta, bandpass_alpha]):   \n        hilbert_node = WrappedPartial(hilbert_transform).getfunc('hilbert_'+band)\n        preproc_tree.add_child_to_node(child, hilbert_node)\n    #     Level 3 functions\n        preproc_tree.add_child_to_node(hilbert_node, WrappedPartial(amplitude).getfunc('amplitude_'+band))\n        preproc_tree.add_child_to_node(hilbert_node, WrappedPartial(phase).getfunc('phase_'+band))\n    preproc_tree.add_child_to_node(average_20_hz,acceleration)\n    preproc_tree.add_child_to_node(average_20_hz, cartesian_to_spherical)\n    preproc_tree.add_child_to_node(average_20_hz, velocity)\n    return preproc_tree\n\npreproc_tree = fog_pipeline(fs=128)\n# Visualize the tree\ngraph = preproc_tree.visualize(size=('16,9'))\ngraph.render(filename='function_tree', directory='.', cleanup=True)\nimg = mpimg.imread('function_tree.png')\nplt.figure(figsize=(12, 5))\nimgplot = plt.imshow(img)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-02-25T11:49:36.413258Z","iopub.execute_input":"2024-02-25T11:49:36.413719Z","iopub.status.idle":"2024-02-25T11:49:37.160511Z","shell.execute_reply.started":"2024-02-25T11:49:36.413687Z","shell.execute_reply":"2024-02-25T11:49:37.159156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. Preprocess and transform all training data","metadata":{}},{"cell_type":"code","source":"session_id = defog_train_sessions[0]\ndf = pd.read_csv(defog_train_dir / '{}.csv'.format(session_id))\ndf","metadata":{"execution":{"iopub.status.busy":"2024-02-25T11:49:40.377365Z","iopub.execute_input":"2024-02-25T11:49:40.377772Z","iopub.status.idle":"2024-02-25T11:49:40.749683Z","shell.execute_reply.started":"2024-02-25T11:49:40.377743Z","shell.execute_reply":"2024-02-25T11:49:40.748311Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# raise\n\ndef prepare_fog_traning(sessions_dir, session_ids, fs, window_size, stride):\n\n    X_all = []\n    y_all = []\n    preproc_tree = fog_pipeline(fs=fs)\n\n    for session_id in tqdm(session_ids):\n        df = pd.read_csv(sessions_dir / '{}.csv'.format(session_id))\n        x_arr = df[['AccV','AccML', 'AccAP']].to_numpy().T\n        y_arr = df[['StartHesitation', 'Turn', 'Walking']].to_numpy().T\n        if 'Valid' in df.columns:\n            mask = np.all(df[['Valid', 'Task']].to_numpy(), axis=1)\n            x_arr = x_arr[:, mask]\n            y_arr = y_arr[:, mask]\n        if np.sum(y_arr) == 0:\n            continue\n        del df\n        gc.collect()\n        x_arr = preproc_tree.evaluate(x_arr)\n        x_arr = x_arr.reshape((x_arr.shape[0]*x_arr.shape[1], -1))\n        X_all.append(x_arr.T)\n        del x_arr\n        gc.collect()\n        y_all.append(y_arr.T)\n        del y_arr\n        gc.collect()\n    \n    X_train = []\n    y_train = []\n    \n    for x_long, y_long in zip(X_all, y_all): \n        for x, y in generate_shorter_sequences(x_long, y_long, window_size, stride):\n            X_train.append(x)\n            y = np.hstack([y, ~np.any(y, axis=1)[..., None]])\n            y_train.append(y)\n        \n    X_train, y_train = np.array(X_train), np.array(y_train)\n        \n    return X_train, y_train\n    ","metadata":{"execution":{"iopub.status.busy":"2024-02-25T11:49:48.389825Z","iopub.execute_input":"2024-02-25T11:49:48.390313Z","iopub.status.idle":"2024-02-25T11:49:48.406530Z","shell.execute_reply.started":"2024-02-25T11:49:48.390279Z","shell.execute_reply":"2024-02-25T11:49:48.403988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_train_tdcsfog, y_train_tdcsfog = prepare_fog_traning(tdcsfog_train_dir, tdcsfog_train_sessions, fs=128, window_size=3406, stride=128*10)\nX_train_defog, y_train_defog = prepare_fog_traning(defog_train_dir, defog_train_sessions, fs=100, window_size=100*60, stride=100*50)","metadata":{"execution":{"iopub.status.busy":"2024-02-25T11:49:52.002659Z","iopub.execute_input":"2024-02-25T11:49:52.004397Z","iopub.status.idle":"2024-02-25T12:01:21.483749Z","shell.execute_reply.started":"2024-02-25T11:49:52.004338Z","shell.execute_reply":"2024-02-25T12:01:21.482276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_points = np.sum(y_train_defog)\nfor i in range(y_train_defog.shape[-1]):\n    ratio = 1/ (y_train_defog[:, :, i].sum() / all_points)\n    print(i, ratio)\n    ","metadata":{"execution":{"iopub.status.busy":"2024-02-24T16:51:32.744258Z","iopub.execute_input":"2024-02-24T16:51:32.744541Z","iopub.status.idle":"2024-02-24T16:51:32.816055Z","shell.execute_reply.started":"2024-02-24T16:51:32.744516Z","shell.execute_reply":"2024-02-24T16:51:32.815137Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. Train the model - tDCS FOG","metadata":{}},{"cell_type":"markdown","source":"## 3.1. Set up","metadata":{}},{"cell_type":"markdown","source":"### 3.1.1. Training config","metadata":{}},{"cell_type":"code","source":"# Define the weights for different classes in the loss function\nloss_weights_tdcsfog = [2, 1, 1, 0.1]\nloss_weights_defog = [500, 7, 41, 1]\n\n# Determine the shape of the input data\ninput_shape_tdcs = (np.shape(X_train_tdcsfog)[-2], np.shape(X_train_tdcsfog)[-1])\ninput_shape_defog = (np.shape(X_train_defog)[-2], np.shape(X_train_defog)[-1])\n\n# Define the optimizer with a specific learning rate\noptimizer_tdcsfog = keras.optimizers.Adam(learning_rate=0.01)\noptimizer_defog = keras.optimizers.Adam(learning_rate=0.01)\n\n# Define a learning rate scheduler to adjust the learning rate based on validation loss\nsheduler_tdcsfog = ReduceLROnPlateau(monitor=\"val_loss\", factor=0.1, patience=10, verbose=1)\nsheduler_defog = ReduceLROnPlateau(monitor=\"val_loss\", factor=0.1, patience=10, verbose=1)\n\n# Define early stopping to prevent overfitting by monitoring validation loss\nearly_stopping_tdcsfog = EarlyStopping(monitor='val_loss', patience=20, restore_best_weights=True, mode='min')\nearly_stopping_defog = EarlyStopping(monitor='val_loss', patience=20, restore_best_weights=True, mode='min')\n\n# Define the neural network model architecture using Keras Sequential API\nmodel_tdcsfog = tf.keras.Sequential([\n    layers.Bidirectional(layers.LSTM(32, return_sequences=True, input_shape=input_shape_tdcs), merge_mode=\"ave\"),  # LSTM layer with 32 units\n    layers.Dense(4, activation=\"softmax\", input_shape=input_shape_tdcs)  # Dense output layer with softmax activation\n])\nmodel_defog = tf.keras.Sequential([\n    layers.Bidirectional(layers.LSTM(32, return_sequences=True, input_shape=input_shape_defog), merge_mode=\"ave\"),  # LSTM layer with 32 units\n    layers.Dense(4, activation=\"softmax\", input_shape=input_shape_defog)  # Dense output layer with softmax activation\n])\n\n# Define a CSV logger to log training history to a file\ncsv_logger_tdcsfog = CSVLogger(\"/kaggle/working/model_history_log_tdcsfog.csv\", append=True)\ncsv_logger_defog = CSVLogger(\"/kaggle/working/model_history_log_defog.csv\", append=True)\n\n# Compile the model with specified optimizer, loss function, and metrics\nmodel_tdcsfog.compile(\n    optimizer=optimizer_tdcsfog,\n    loss='categorical_crossentropy',  # Categorical crossentropy loss function\n    loss_weights=loss_weights_tdcsfog,  # Apply class weights to the loss function\n    weighted_metrics=['categorical_accuracy'],# Use weighted accuracy as a metric\n    metrics=[tf.keras.metrics.Precision(class_id=0),\n            tf.keras.metrics.Precision(class_id=1),\n            tf.keras.metrics.Precision(class_id=2),\n            tf.keras.metrics.Recall(class_id=0),\n            tf.keras.metrics.Recall(class_id=1),\n            tf.keras.metrics.Recall(class_id=2)\n            ]\n)\nmodel_defog.compile(\n    optimizer=optimizer_defog,\n    loss='categorical_crossentropy',  # Categorical crossentropy loss function\n    loss_weights=loss_weights_defog,  # Apply class weights to the loss function\n    weighted_metrics=['categorical_accuracy'],# Use weighted accuracy as a metric\n    metrics=[tf.keras.metrics.Precision(class_id=0),\n            tf.keras.metrics.Precision(class_id=1),\n            tf.keras.metrics.Precision(class_id=2),\n            tf.keras.metrics.Recall(class_id=0),\n            tf.keras.metrics.Recall(class_id=1),\n            tf.keras.metrics.Recall(class_id=2)\n            ]\n)\n\n","metadata":{"execution":{"iopub.status.busy":"2024-02-24T16:51:32.818786Z","iopub.execute_input":"2024-02-24T16:51:32.819066Z","iopub.status.idle":"2024-02-24T16:51:34.322519Z","shell.execute_reply.started":"2024-02-24T16:51:32.819042Z","shell.execute_reply":"2024-02-24T16:51:34.321721Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 3.1.2. Data preparation","metadata":{}},{"cell_type":"code","source":"# Set seed for reproducibility\nSEED = 42\nnp.random.seed(SEED)\n\n# Define the validation split ratio\nval_split = 0.1\n\n# Calculate the number of samples for validation\nval_samples_tdcsfog = int(len(X_train_tdcsfog) * val_split)\nval_samples_defog = int(len(X_train_defog) * val_split)\n\n# Randomly select samples for validation from the entire dataset\nval_id_tdcsfog = np.random.choice(len(X_train_tdcsfog), val_samples_tdcsfog, replace=False)\nval_id_defog = np.random.choice(len(X_train_defog), val_samples_defog, replace=False)\n\n# Create training and validation datasets using indices\ntrain_id_tdcsfog = np.setdiff1d(np.arange(len(X_train_tdcsfog)), val_id_tdcsfog)\ntrain_dataset_tdcsfog = tf.data.Dataset.from_tensor_slices((X_train_tdcsfog[train_id_tdcsfog], y_train_tdcsfog[train_id_tdcsfog]))\nval_dataset_tdcsfog = tf.data.Dataset.from_tensor_slices((X_train_tdcsfog[val_id_tdcsfog], y_train_tdcsfog[val_id_tdcsfog]))\n\ntrain_id_defog = np.setdiff1d(np.arange(len(X_train_defog)), val_id_defog)\ntrain_dataset_defog = tf.data.Dataset.from_tensor_slices((X_train_defog[train_id_defog], y_train_defog[train_id_defog]))\nval_dataset_defog = tf.data.Dataset.from_tensor_slices((X_train_defog[val_id_defog], y_train_defog[val_id_defog]))\n\n# Define the batch size\nbatch_size = 32\n\n# Shuffle and batch the training dataset\ntrain_dataset_tdcsfog = train_dataset_tdcsfog.shuffle(buffer_size=len(train_id_tdcsfog),\n                                      seed=SEED).batch(batch_size)\ntrain_dataset_defog = train_dataset_defog.shuffle(buffer_size=len(train_id_defog),\n                                      seed=SEED).batch(batch_size)\n\n# Batch the validation dataset\nval_dataset_tdcsfog = val_dataset_tdcsfog.batch(batch_size)\nval_dataset_defog = val_dataset_defog.batch(batch_size)","metadata":{"execution":{"iopub.status.busy":"2024-02-24T16:51:34.323664Z","iopub.execute_input":"2024-02-24T16:51:34.324014Z","iopub.status.idle":"2024-02-24T16:51:43.840681Z","shell.execute_reply.started":"2024-02-24T16:51:34.323979Z","shell.execute_reply":"2024-02-24T16:51:43.83989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.2. Training routine","metadata":{}},{"cell_type":"code","source":"history_tdcsfog = model_tdcsfog.fit(\n    train_dataset_tdcsfog,\n    epochs=200,\n    validation_data=val_dataset_tdcsfog,\n    verbose=1,\n    callbacks=[early_stopping_tdcsfog, sheduler_tdcsfog, csv_logger_tdcsfog],\n#     use_multiprocessing=True\n)\n\n# Store training history as a dataframe\nhistory_df_tdcsfog = pd.DataFrame(history_tdcsfog.history)","metadata":{"execution":{"iopub.status.busy":"2024-02-24T16:51:43.841775Z","iopub.execute_input":"2024-02-24T16:51:43.842066Z","iopub.status.idle":"2024-02-24T16:52:23.925082Z","shell.execute_reply.started":"2024-02-24T16:51:43.842041Z","shell.execute_reply":"2024-02-24T16:52:23.924345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history_defog = model_defog.fit(\n    train_dataset_defog,\n    epochs=200,\n    validation_data=val_dataset_defog,\n    verbose=1,\n    callbacks=[early_stopping_defog, sheduler_defog, csv_logger_defog],\n#     use_multiprocessing=True\n)\n\n# Store training history as a dataframe\nhistory_df_defog = pd.DataFrame(history_defog.history)","metadata":{"execution":{"iopub.status.busy":"2024-02-24T16:52:23.926195Z","iopub.execute_input":"2024-02-24T16:52:23.926753Z","iopub.status.idle":"2024-02-24T16:52:42.609649Z","shell.execute_reply.started":"2024-02-24T16:52:23.926722Z","shell.execute_reply":"2024-02-24T16:52:42.608717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. Submission","metadata":{}},{"cell_type":"markdown","source":"## 4.1. Preprocess the test set","metadata":{}},{"cell_type":"code","source":"del train_dataset_tdcsfog, train_dataset_defog, val_dataset_tdcsfog, val_dataset_defog, X_train_tdcsfog, y_train_tdcsfog, X_train_defog, y_train_defog\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-02-24T17:30:28.146844Z","iopub.execute_input":"2024-02-24T17:30:28.147538Z","iopub.status.idle":"2024-02-24T17:30:28.492437Z","shell.execute_reply.started":"2024-02-24T17:30:28.147508Z","shell.execute_reply":"2024-02-24T17:30:28.490984Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def prepare_fog_testing(sessions_dir, session_ids, fs, window_size, stride):\n    X_all = []\n    X_ids = []\n    preproc_tree = fog_pipeline(fs=fs)\n\n    for session_id in tqdm(session_ids):\n        df = pd.read_csv(sessions_dir / '{}.csv'.format(session_id))\n        id_time = [f'{session_id}_{t}' for t in range(len(df))]\n        X_ids.append(np.array(id_time)[..., None])\n        x_arr = df[['AccV','AccML', 'AccAP']].to_numpy().T\n        if 'Valid' in df.columns:\n            mask = np.all(df[['Valid', 'Task']].to_numpy(), axis=1)\n            x_arr = x_arr[:, mask]\n        del df\n        gc.collect()\n        x_arr = preproc_tree.evaluate(x_arr)\n        x_arr = x_arr.reshape((x_arr.shape[0]*x_arr.shape[1], -1))\n        X_all.append(x_arr.T)\n        del x_arr\n        gc.collect()\n    \n    X_test = []\n    X_id_time = []\n    \n    for x_long, id_long in zip(X_all, X_ids): \n        for x, i in generate_shorter_sequences(x_long, id_long, window_size, stride):\n            X_test.append(x)\n            X_id_time.append(i)\n            \n        \n    X_test, X_id_time = np.array(X_test), np.array(X_id_time)\n    if X_test.ndim == 2:\n        X_test = X_test[None, ...]\n        \n    return X_test, X_id_time\n\n\ndef create_template(ids, example=sample_submission):\n    template = pd.DataFrame(np.zeros((len(ids), 4)), columns=sample_submission.columns)\n    template['Id'] = ids\n    return template\n\n\ndef populate_submission(model, test_dataset, template):\n    for batch, id_batch in test_dataset:\n        pred = model.predict(batch, verbose=1, use_multiprocessing=True)\n        for seqi, id_single in enumerate(id_batch):\n            first_id = id_single[0]\n            first_id = tf.compat.as_str_any(first_id.numpy()[0])\n            last_id = id_single[-1]\n            last_id = tf.compat.as_str_any(last_id.numpy()[0])\n            \n            first_loc = template.loc[template['Id']==first_id].index.tolist()[0]\n            last_loc = template.loc[template['Id']==last_id].index.tolist()[0]\n              \n            template.iloc[first_loc:last_loc+1, 1:] = pred[seqi, :, :-1]\n            \n#             for row, t in tqdm(enumerate(i)):\n#                 tval = tf.compat.as_str_any(t.numpy()[0])\n#                 if template[template['Id']==tval].iloc[:, 1:].sum(axis=1).iloc[0] != 0:\n#                     template[template['Id']==tval] = np.hstack((np.array(tval), pred[seqi, row, :-1].ravel()))\n#                 if tval not in existing_id:\n#                     submission_csv.loc[len(submission_csv)] = np.hstack((np.array(tval), pred[seqi, row, :-1].ravel()))\n    return template","metadata":{"execution":{"iopub.status.busy":"2024-02-24T17:24:47.79867Z","iopub.execute_input":"2024-02-24T17:24:47.799033Z","iopub.status.idle":"2024-02-24T17:24:47.816395Z","shell.execute_reply.started":"2024-02-24T17:24:47.799003Z","shell.execute_reply":"2024-02-24T17:24:47.815495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ids = []\nfor session in tqdm(tdcsfog_test_sessions):\n    df = pd.read_csv(tdcsfog_test_dir / '{}.csv'.format(session))\n    this_ids = [str(session)+'_'+str(t) for t in range(len(df))]\n    ids = ids + this_ids\n        \nfor session in tqdm(defog_test_sessions):\n    df = pd.read_csv(defog_test_dir / '{}.csv'.format(session))\n    this_ids = [str(session)+'_'+str(t) for t in range(len(df))]\n    ids = ids + this_ids\n        \ntemplate = create_template(ids)","metadata":{"execution":{"iopub.status.busy":"2024-02-24T17:10:51.21359Z","iopub.execute_input":"2024-02-24T17:10:51.214303Z","iopub.status.idle":"2024-02-24T17:10:51.644158Z","shell.execute_reply.started":"2024-02-24T17:10:51.21426Z","shell.execute_reply":"2024-02-24T17:10:51.643191Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = None\nfor session in tdcsfog_test_sessions:\n    X_test_tdcsfog, X_test_tdcsfog_idt = prepare_fog_testing(tdcsfog_test_dir, [session], fs=128, window_size=3406, stride=128*24)\n    X_test_dataset_tdcsfog = tf.data.Dataset.from_tensor_slices((X_test_tdcsfog, X_test_tdcsfog_idt)).batch(32)\n    test_df_tdcsfog = populate_submission(model_tdcsfog, X_test_dataset_tdcsfog, template)\n    del test_df_tdcsfog\n    gc.collect()\n\nfor session in defog_test_sessions:\n    X_test_defog, X_test_defog_idt = prepare_fog_testing(defog_test_dir, [session], fs=100, window_size=100*60, stride=100*58)\n    X_test_dataset_defog = tf.data.Dataset.from_tensor_slices((X_test_defog, X_test_defog_idt)).batch(32)\n    test_df_defog = populate_submission(model_defog, X_test_dataset_defog, template)\n    del test_df_defog\n    gc.collect()\n    \ntemplate.to_csv('submission.csv', index=False)\ntemplate","metadata":{"execution":{"iopub.status.busy":"2024-02-24T17:25:04.994858Z","iopub.execute_input":"2024-02-24T17:25:04.995796Z","iopub.status.idle":"2024-02-24T17:25:19.850701Z","shell.execute_reply.started":"2024-02-24T17:25:04.995759Z","shell.execute_reply":"2024-02-24T17:25:19.849738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# model.save('/kaggle/working/exp4_lstm/')\n# import shutil\n# shutil.make_archive('exp3_lstm', 'zip', '/kaggle/working/exp4_lstm/')","metadata":{"execution":{"iopub.status.busy":"2024-02-20T10:42:36.254387Z","iopub.execute_input":"2024-02-20T10:42:36.255304Z","iopub.status.idle":"2024-02-20T10:42:39.529065Z","shell.execute_reply.started":"2024-02-20T10:42:36.255269Z","shell.execute_reply":"2024-02-20T10:42:39.528237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Scratchpad","metadata":{}},{"cell_type":"markdown","source":"# Pure data","metadata":{}},{"cell_type":"code","source":"# sample_sessions = '76c7edf878'\n# sample_df = pd.read_csv(tdcsfog_train_dir / '{}.csv'.format(sample_sessions))\n# for column in ['AccV','AccML', 'AccAP']:\n#     plt.figure(figsize=(10, 1))\n#     for color, event in zip(['red', 'blue', 'green'], \n#                             ['StartHesitation', 'Turn', 'Walking']):\n#         mask_event = np.argwhere(sample_df[event].to_numpy().astype(bool)).ravel()\n#         timeseries_event = sample_df[column]\n#         sum_acc = timeseries_event.sum()\n#         timepoints = np.arange(len(timeseries_event))\n#         mean_acc = timepoints[mask_event].mean()\n#         plt.plot(timepoints[mask_event] / 128, timeseries_event[mask_event], color=color, label=event+', mean acc = {:.4f}'.format(mean_acc))\n#         plt.title(column)\n#         plt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-12-19T18:27:42.20375Z","iopub.execute_input":"2023-12-19T18:27:42.204819Z","iopub.status.idle":"2023-12-19T18:27:44.098606Z","shell.execute_reply.started":"2023-12-19T18:27:42.204779Z","shell.execute_reply":"2023-12-19T18:27:44.097666Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sample_sessions = '76c7edf878'\n# sample_df = pd.read_csv(tdcsfog_train_dir / '{}.csv'.format(sample_sessions))\n# window_size = 128 * 4\n# for column in ['AccV','AccML', 'AccAP']:\n#     plt.figure(figsize=(10, 1))\n#     for color, event in zip(['red', 'blue', 'green'], \n#                             ['StartHesitation', 'Turn', 'Walking']):\n#         mask_event = np.argwhere(sample_df[event].to_numpy().astype(bool)).ravel()\n#         timeseries_event = sample_df[column]\n#         sum_acc = timeseries_event.sum()\n#         timepoints = np.arange(len(timeseries_event))\n#         x =  timepoints[mask_event]\n#         y = timeseries_event[mask_event].rolling(window=window_size).mean()\n#         plt.plot(x / 128, y, color=color, label=event)\n#         plt.title(column)\n#         plt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-12-21T10:38:12.323084Z","iopub.execute_input":"2023-12-21T10:38:12.323547Z","iopub.status.idle":"2023-12-21T10:38:14.490906Z","shell.execute_reply.started":"2023-12-21T10:38:12.323481Z","shell.execute_reply":"2023-12-21T10:38:14.489824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# ","metadata":{}},{"cell_type":"markdown","source":"# Spherical","metadata":{}},{"cell_type":"code","source":"# sf = 128\n# start, stop = 5, 500 # seconds\n# start_bin, stop_bin = start * sf, stop * sf\n# sample_subject = '76c7edf878'\n# sample_df = pd.read_csv(tdcsfog_train_dir / '{}.csv'.format(sample_subject))\n\n# r_list = []\n# theta_list = []\n# phi_list = []\n\n# # compute the radial distance, azimuth angle, polar angle\n# for vector in sample_df[['AccV','AccML', 'AccAP']].to_numpy():\n#     r, theta, phi = cartesian_to_spherical(*vector)\n#     r_list.append(r)\n#     theta_list.append(theta)\n#     phi_list.append(phi)\n\n# for name, unit, feature in zip(['vector norm', 'polar angle', 'azimuth angle'],\n#                                ['G', 'radian', 'radian'],\n#                          [np.array(r_list), np.array(theta_list), np.array(phi_list)]):\n#     plt.figure(figsize=(10, 1))\n#     for color, event in zip(['red', 'blue', 'green'], \n#                             ['StartHesitation', 'Turn', 'Walking']):\n#         mask_event = sample_df[event].to_numpy().astype(bool)\n#         timeseries_event = feature\n#         timepoints = np.arange(len(timeseries_event))\n#         epoch_mask = (start_bin <= timepoints) & (timepoints <= stop_bin)\n#         mask = epoch_mask & mask_event\n#         plt.plot(timepoints[mask] / 128, timeseries_event[mask], color=color, label=event)\n#         plt.title(name)\n#         plt.ylabel(unit)\n#         plt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-12-19T18:16:26.806795Z","iopub.execute_input":"2023-12-19T18:16:26.807296Z","iopub.status.idle":"2023-12-19T18:16:29.533064Z","shell.execute_reply.started":"2023-12-19T18:16:26.807257Z","shell.execute_reply":"2023-12-19T18:16:29.531977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Correlation between cartesian and spherical accelerations","metadata":{}},{"cell_type":"code","source":"# accelerations = sample_df[['AccV','AccML', 'AccAP']].to_numpy()\n# for ft in [r_list, theta_list, phi_list]:\n#     ft = np.array(ft)[..., None]\n#     accelerations = np.hstack([accelerations, ft])\n    \n# print(np.corrcoef(accelerations.T))\n# # accelerations.shape\n    ","metadata":{"execution":{"iopub.status.busy":"2023-12-19T17:33:05.734245Z","iopub.execute_input":"2023-12-19T17:33:05.734744Z","iopub.status.idle":"2023-12-19T17:33:05.762585Z","shell.execute_reply.started":"2023-12-19T17:33:05.734715Z","shell.execute_reply":"2023-12-19T17:33:05.761649Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Speed","metadata":{}},{"cell_type":"code","source":"# sf = 128\n# start, stop = 0, 500 # seconds\n# start_bin, stop_bin = start * sf, stop * sf\n\n# sample_sessions = '76c7edf878'\n# sample_df = pd.read_csv(tdcsfog_train_dir / '{}.csv'.format(sample_sessions))\n# for column in ['AccV','AccML', 'AccAP']:\n#     plt.figure(figsize=(10, 1))\n#     for color, event in zip(['red', 'blue', 'green'], \n#                             ['StartHesitation', 'Turn', 'Walking']):\n        \n# #         mask_event = sample_df[event].to_numpy().astype(bool)\n# #         timeseries_event = feature\n# #         timepoints = np.arange(len(timeseries_event))\n# #         epoch_mask = (start_bin <= timepoints) & (timepoints <= stop_bin)\n# #         mask = epoch_mask & mask_event\n        \n#         mask_event = sample_df[event].to_numpy().astype(bool)\n#         timeseries_event = np.cumsum(sample_df[column])\n#         timepoints = np.arange(len(timeseries_event))\n#         epoch_mask = (start_bin <= timepoints) & (timepoints <= stop_bin)\n#         mask = epoch_mask & mask_event\n#         plt.plot(timepoints[mask] / 128, timeseries_event[mask], color=color, label=event)\n#         plt.title(column)\n#         plt.ylim(-500000, 500000)\n#         plt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-12-19T18:15:49.911664Z","iopub.execute_input":"2023-12-19T18:15:49.912405Z","iopub.status.idle":"2023-12-19T18:15:51.61395Z","shell.execute_reply.started":"2023-12-19T18:15:49.912365Z","shell.execute_reply":"2023-12-19T18:15:51.612836Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Jerk","metadata":{}},{"cell_type":"code","source":"# sf = 128\n# start, stop = 0, 500 # seconds\n# start_bin, stop_bin = start * sf, stop * sf\n\n# sample_sessions = '76c7edf878'\n# sample_df = pd.read_csv(tdcsfog_train_dir / '{}.csv'.format(sample_sessions))\n# for column in ['AccV','AccML', 'AccAP']:\n#     plt.figure(figsize=(10, 1))\n#     for color, event in zip(['red', 'blue', 'green'], \n#                             ['StartHesitation', 'Turn', 'Walking']):\n        \n# #         mask_event = sample_df[event].to_numpy().astype(bool)\n# #         timeseries_event = feature\n# #         timepoints = np.arange(len(timeseries_event))\n# #         epoch_mask = (start_bin <= timepoints) & (timepoints <= stop_bin)\n# #         mask = epoch_mask & mask_event\n        \n#         mask_event = sample_df[event].to_numpy().astype(bool)[1:]\n#         timeseries_event = np.diff(sample_df[column])\n#         timepoints = np.arange(len(timeseries_event))\n#         epoch_mask = (start_bin <= timepoints) & (timepoints <= stop_bin)\n#         mask = epoch_mask & mask_event\n#         plt.plot(timepoints[mask] / 128, timeseries_event[mask], color=color, label=event)\n#         plt.title(column)\n# #         plt.ylim(-500000, 500000)\n#         plt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-12-19T18:14:16.737523Z","iopub.execute_input":"2023-12-19T18:14:16.737999Z","iopub.status.idle":"2023-12-19T18:14:18.719811Z","shell.execute_reply.started":"2023-12-19T18:14:16.737957Z","shell.execute_reply":"2023-12-19T18:14:18.718741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Spectrogram","metadata":{}},{"cell_type":"code","source":"# sf = 128\n# start, stop = 0, 500 # seconds\n# start_bin, stop_bin = start * sf, stop * sf\n\n# sample_sessions = '76c7edf878'\n# sample_df = pd.read_csv(tdcsfog_train_dir / '{}.csv'.format(sample_sessions))\n# for column in ['AccV','AccML', 'AccAP']:\n#     plt.figure(figsize=(10, 1))\n#     for color, event in zip(['red', 'blue', 'green'], \n#                             ['StartHesitation', 'Turn', 'Walking']):\n        \n# #         mask_event = sample_df[event].to_numpy().astype(bool)\n# #         timeseries_event = feature\n# #         timepoints = np.arange(len(timeseries_event))\n# #         epoch_mask = (start_bin <= timepoints) & (timepoints <= stop_bin)\n# #         mask = epoch_mask & mask_event\n        \n#         mask_event = sample_df[event].to_numpy().astype(bool)[1:]\n#         timeseries_event = np.diff(sample_df[column])\n#         timepoints = np.arange(len(timeseries_event))\n#         epoch_mask = (start_bin <= timepoints) & (timepoints <= stop_bin)\n#         mask = epoch_mask & mask_event\n#         plot_spectrogram(timeseries_event[mask], fs=100, maxfreq=50, minfreq=0)\n# #         plt.plot(timepoints[mask] / 128, timeseries_event[mask], color=color, label=event)\n#         plt.title(column + ' ' + event)\n# # #         plt.ylim(-500000, 500000)\n# #         plt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-12-19T18:59:58.78109Z","iopub.execute_input":"2023-12-19T18:59:58.781469Z","iopub.status.idle":"2023-12-19T19:00:01.444837Z","shell.execute_reply.started":"2023-12-19T18:59:58.78144Z","shell.execute_reply":"2023-12-19T19:00:01.443663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Hilbert","metadata":{}},{"cell_type":"code","source":"# # Bandpass filter parameters\n# lowcut = 1  # Lower cutoff frequency\n# highcut = 5  # Higher cutoff frequency\n# fs = 100  # Sampling frequency\n\n\n# sf = 128\n# start, stop = 0, 500 # seconds\n# start_bin, stop_bin = start * sf, stop * sf\n\n# sample_sessions = '76c7edf878'\n# sample_df = pd.read_csv(tdcsfog_train_dir / '{}.csv'.format(sample_sessions))\n# for column in ['AccV','AccML', 'AccAP']:\n#     plt.figure(figsize=(10, 1))\n#     for color, event in zip(['red', 'blue', 'green'], \n#                             ['StartHesitation', 'Turn', 'Walking']):\n        \n#         mask_event = sample_df[event].to_numpy().astype(bool)[1:]\n#         timeseries_event = np.diff(sample_df[column])\n#         timepoints = np.arange(len(timeseries_event))\n#         epoch_mask = (start_bin <= timepoints) & (timepoints <= stop_bin)\n#         mask = epoch_mask & mask_event\n#         x = timepoints[mask] / 128\n#         y = timeseries_event[mask]\n#         fig, ax = plt.subplots(2, 1, sharex=True)\n#         # Apply bandpass filter\n#         filtered_data = butter_bandpass_filter(y, lowcut, highcut, fs)\n\n#         # Compute Hilbert transform\n#         analytic_signal = hilbert(filtered_data)\n#         amplitude_envelope = np.abs(analytic_signal)\n#         instantaneous_phase = np.angle(analytic_signal)\n        \n#         ax[0].plot(instantaneous_phase, color=color)\n#         ax[0].set_title('Phase, '+event)\n#         ax[1].plot(amplitude_envelope, color=color)\n#         ax[1].set_title('Amplitude, '+event)\n        \n#         plt.title(column)\n# #         plt.ylim(-500000, 500000)\n#         plt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-12-19T19:15:33.661369Z","iopub.execute_input":"2023-12-19T19:15:33.661777Z","iopub.status.idle":"2023-12-19T19:15:39.239895Z","shell.execute_reply.started":"2023-12-19T19:15:33.661745Z","shell.execute_reply":"2023-12-19T19:15:39.238942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}