.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/initial_contacts/_03_hklee_algo.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_initial_contacts__03_hklee_algo.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_initial_contacts__03_hklee_algo.py:


HKLee algo
==========

This example shows how to use the improved HKLee algorithm and some examples on how the results compare to the original
matlab implementation.

.. GENERATED FROM PYTHON SOURCE LINES 9-15

.. code-block:: default


    import pandas as pd
    from matplotlib import pyplot as plt
    from mobgap.data import LabExampleDataset
    from mobgap.initial_contacts._hklee_algo_improved import IcdHKLeeImproved








.. GENERATED FROM PYTHON SOURCE LINES 16-25

Loading data
------------
.. note :: More infos about data loading can be found in the :ref:`data loading example <data_loading_example>`.

We load example data from the lab dataset together with the INDIP reference system.
We will use the INDIP output for initial contacts ("ic") as ground truth.

Like most algorithms, the GsdIluz requires the data to be in body frame coordinates.
As we know the sensor was well aligned, we can just use ``to_body_frame`` to transform the data.

.. GENERATED FROM PYTHON SOURCE LINES 25-27

.. code-block:: default

    from mobgap.utils.conversions import to_body_frame








.. GENERATED FROM PYTHON SOURCE LINES 28-33

Loading some data
-----------------
We load example data from the lab dataset together with the INDIP reference system.
Like most algorithms, the algorithm requires the data to be in body frame coordinates.
As we know the sensor was well aligned, we can just use ``to_body_frame`` to transform the data.

.. GENERATED FROM PYTHON SOURCE LINES 33-48

.. code-block:: default


    example_data = LabExampleDataset(
        reference_system="INDIP", reference_para_level="wb"
    )

    single_test = example_data.get_subset(
        cohort="HA", participant_id="001", test="Test11", trial="Trial1"
    )
    imu_data = to_body_frame(single_test.data_ss)
    reference_wbs = single_test.reference_parameters_.wb_list

    sampling_rate_hz = single_test.sampling_rate_hz
    ref_ics = single_test.reference_parameters_.ic_list

    reference_wbs





.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>start</th>
          <th>end</th>
          <th>n_strides</th>
          <th>duration_s</th>
          <th>length_m</th>
          <th>avg_walking_speed_mps</th>
          <th>avg_cadence_spm</th>
          <th>avg_stride_length_m</th>
          <th>termination_reason</th>
        </tr>
        <tr>
          <th>wb_id</th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>632</td>
          <td>988</td>
          <td>5</td>
          <td>3.55</td>
          <td>3.428989</td>
          <td>0.975373</td>
          <td>104.069084</td>
          <td>1.124391</td>
          <td>Pause</td>
        </tr>
        <tr>
          <th>1</th>
          <td>2864</td>
          <td>3325</td>
          <td>4</td>
          <td>4.60</td>
          <td>1.452572</td>
          <td>0.411857</td>
          <td>81.296475</td>
          <td>0.581029</td>
          <td>Pause</td>
        </tr>
        <tr>
          <th>2</th>
          <td>3853</td>
          <td>5085</td>
          <td>16</td>
          <td>12.31</td>
          <td>7.044042</td>
          <td>0.617801</td>
          <td>89.246331</td>
          <td>0.838960</td>
          <td>Pause</td>
        </tr>
        <tr>
          <th>3</th>
          <td>7641</td>
          <td>8621</td>
          <td>12</td>
          <td>9.79</td>
          <td>4.396574</td>
          <td>0.510108</td>
          <td>94.370318</td>
          <td>0.645176</td>
          <td>Pause</td>
        </tr>
        <tr>
          <th>4</th>
          <td>9451</td>
          <td>9932</td>
          <td>6</td>
          <td>4.80</td>
          <td>3.545277</td>
          <td>0.755728</td>
          <td>88.778698</td>
          <td>1.021695</td>
          <td>Pause</td>
        </tr>
        <tr>
          <th>5</th>
          <td>11989</td>
          <td>12517</td>
          <td>6</td>
          <td>5.27</td>
          <td>3.514735</td>
          <td>0.880632</td>
          <td>95.832693</td>
          <td>1.021576</td>
          <td>Pause</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 49-53

Applying the algorithm
----------------------
Below we apply the shin algorithm to a lab trial.
We will use the `GsIterator` to iterate over the gait sequences and apply the algorithm to each wb.

.. GENERATED FROM PYTHON SOURCE LINES 53-67

.. code-block:: default

    from mobgap.pipeline import GsIterator

    iterator = GsIterator()

    for (gs, data), result in iterator.iterate(imu_data, reference_wbs):
        result.ic_list = (
            IcdHKLeeImproved()
            .detect(data, sampling_rate_hz=sampling_rate_hz)
            .ic_list_
        )

    detected_ics = iterator.results_.ic_list

    detected_ics





.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th></th>
          <th>ic</th>
        </tr>
        <tr>
          <th>wb_id</th>
          <th>step_id</th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th rowspan="6" valign="top">0</th>
          <th>0</th>
          <td>649</td>
        </tr>
        <tr>
          <th>1</th>
          <td>712</td>
        </tr>
        <tr>
          <th>2</th>
          <td>772</td>
        </tr>
        <tr>
          <th>3</th>
          <td>829</td>
        </tr>
        <tr>
          <th>4</th>
          <td>885</td>
        </tr>
        <tr>
          <th>5</th>
          <td>941</td>
        </tr>
        <tr>
          <th rowspan="7" valign="top">1</th>
          <th>0</th>
          <td>2879</td>
        </tr>
        <tr>
          <th>1</th>
          <td>2943</td>
        </tr>
        <tr>
          <th>2</th>
          <td>3006</td>
        </tr>
        <tr>
          <th>3</th>
          <td>3075</td>
        </tr>
        <tr>
          <th>4</th>
          <td>3127</td>
        </tr>
        <tr>
          <th>5</th>
          <td>3204</td>
        </tr>
        <tr>
          <th>6</th>
          <td>3276</td>
        </tr>
        <tr>
          <th rowspan="17" valign="top">2</th>
          <th>0</th>
          <td>3866</td>
        </tr>
        <tr>
          <th>1</th>
          <td>3927</td>
        </tr>
        <tr>
          <th>2</th>
          <td>3994</td>
        </tr>
        <tr>
          <th>3</th>
          <td>4093</td>
        </tr>
        <tr>
          <th>4</th>
          <td>4156</td>
        </tr>
        <tr>
          <th>5</th>
          <td>4225</td>
        </tr>
        <tr>
          <th>6</th>
          <td>4292</td>
        </tr>
        <tr>
          <th>7</th>
          <td>4365</td>
        </tr>
        <tr>
          <th>8</th>
          <td>4438</td>
        </tr>
        <tr>
          <th>9</th>
          <td>4504</td>
        </tr>
        <tr>
          <th>10</th>
          <td>4573</td>
        </tr>
        <tr>
          <th>11</th>
          <td>4651</td>
        </tr>
        <tr>
          <th>12</th>
          <td>4706</td>
        </tr>
        <tr>
          <th>13</th>
          <td>4771</td>
        </tr>
        <tr>
          <th>14</th>
          <td>4827</td>
        </tr>
        <tr>
          <th>15</th>
          <td>4879</td>
        </tr>
        <tr>
          <th>16</th>
          <td>4984</td>
        </tr>
        <tr>
          <th rowspan="15" valign="top">3</th>
          <th>0</th>
          <td>7654</td>
        </tr>
        <tr>
          <th>1</th>
          <td>7748</td>
        </tr>
        <tr>
          <th>2</th>
          <td>7797</td>
        </tr>
        <tr>
          <th>3</th>
          <td>7858</td>
        </tr>
        <tr>
          <th>4</th>
          <td>7913</td>
        </tr>
        <tr>
          <th>5</th>
          <td>8003</td>
        </tr>
        <tr>
          <th>6</th>
          <td>8059</td>
        </tr>
        <tr>
          <th>7</th>
          <td>8111</td>
        </tr>
        <tr>
          <th>8</th>
          <td>8181</td>
        </tr>
        <tr>
          <th>9</th>
          <td>8241</td>
        </tr>
        <tr>
          <th>10</th>
          <td>8293</td>
        </tr>
        <tr>
          <th>11</th>
          <td>8352</td>
        </tr>
        <tr>
          <th>12</th>
          <td>8414</td>
        </tr>
        <tr>
          <th>13</th>
          <td>8480</td>
        </tr>
        <tr>
          <th>14</th>
          <td>8557</td>
        </tr>
        <tr>
          <th rowspan="7" valign="top">4</th>
          <th>0</th>
          <td>9481</td>
        </tr>
        <tr>
          <th>1</th>
          <td>9546</td>
        </tr>
        <tr>
          <th>2</th>
          <td>9609</td>
        </tr>
        <tr>
          <th>3</th>
          <td>9674</td>
        </tr>
        <tr>
          <th>4</th>
          <td>9741</td>
        </tr>
        <tr>
          <th>5</th>
          <td>9809</td>
        </tr>
        <tr>
          <th>6</th>
          <td>9877</td>
        </tr>
        <tr>
          <th rowspan="8" valign="top">5</th>
          <th>0</th>
          <td>12002</td>
        </tr>
        <tr>
          <th>1</th>
          <td>12057</td>
        </tr>
        <tr>
          <th>2</th>
          <td>12114</td>
        </tr>
        <tr>
          <th>3</th>
          <td>12170</td>
        </tr>
        <tr>
          <th>4</th>
          <td>12226</td>
        </tr>
        <tr>
          <th>5</th>
          <td>12289</td>
        </tr>
        <tr>
          <th>6</th>
          <td>12379</td>
        </tr>
        <tr>
          <th>7</th>
          <td>12474</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 68-71

Matlab Outputs
--------------
To check if the algorithm was implemented correctly, we compare the results to the matlab implementation.

.. GENERATED FROM PYTHON SOURCE LINES 71-101

.. code-block:: default

    import json

    from mobgap import PACKAGE_ROOT


    def load_matlab_output(datapoint):
        p = datapoint.group_label
        with (
            PACKAGE_ROOT.parent
            / f"example_data/original_results/icd_hklee_improved/lab/{p.cohort}/{p.participant_id}/SD_Output.json"
        ).open() as f:
            original_results = json.load(f)["SD_Output"][p.time_measure][p.test][
                p.trial
            ]["SU"]["LowerBack"]["SD"]

        if not isinstance(original_results, list):
            original_results = [original_results]

        ics = {}
        for i, gs in enumerate(original_results, start=1):
            ics[i] = pd.DataFrame({"ic": gs["IC"]}).rename_axis(index="step_id")

        return (
            pd.concat(ics, names=["wb_id", ics[1].index.name])
            * datapoint.sampling_rate_hz
        ).astype("int64")


    detected_ics_matlab = load_matlab_output(single_test)
    detected_ics_matlab





.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th></th>
          <th>ic</th>
        </tr>
        <tr>
          <th>wb_id</th>
          <th>step_id</th>
          <th></th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th rowspan="5" valign="top">1</th>
          <th>0</th>
          <td>652</td>
        </tr>
        <tr>
          <th>1</th>
          <td>717</td>
        </tr>
        <tr>
          <th>2</th>
          <td>777</td>
        </tr>
        <tr>
          <th>3</th>
          <td>834</td>
        </tr>
        <tr>
          <th>4</th>
          <td>892</td>
        </tr>
        <tr>
          <th>...</th>
          <th>...</th>
          <td>...</td>
        </tr>
        <tr>
          <th rowspan="5" valign="top">6</th>
          <th>4</th>
          <td>12231</td>
        </tr>
        <tr>
          <th>5</th>
          <td>12294</td>
        </tr>
        <tr>
          <th>6</th>
          <td>12384</td>
        </tr>
        <tr>
          <th>7</th>
          <td>12479</td>
        </tr>
        <tr>
          <th>8</th>
          <td>12514</td>
        </tr>
      </tbody>
    </table>
    <p>64 rows × 1 columns</p>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 102-119

Plotting the results
--------------------
With that we can compare the python, matlab and ground truth results.
We zoom in into one of the gait sequences to better see the output.

We can make a couple of main observations:

1. The python version finds (for the most part) the same ICs as the matlab version, but detects them slightly earlies
   (around 2-8 samples/20-80 ms).
   This is likely due to some differences in the downsampling process.
2. Matlab and Python sometimes differ in the detection in the first and the last IC, if they are close to the border
   of the GS.
   This is because the morphological filter in the Python and the Matlab treat the border differently.
3. Compared to the ground truth reference, both versions detect the IC too late most of the time.
4. Both algorithms can not detect the first IC of the gait sequence.
   However, this is expected, as per definition, this first IC marks the start of the WB in the reference system.
   Hence, there are no samples before that point the algorithm can use to detect the IC.

.. GENERATED FROM PYTHON SOURCE LINES 119-141

.. code-block:: default


    imu_data.reset_index(drop=True).plot(y="acc_is")

    plt.plot(
        ref_ics["ic"], imu_data["acc_is"].iloc[ref_ics["ic"]], "o", label="ref"
    )
    plt.plot(
        detected_ics["ic"],
        imu_data["acc_is"].iloc[detected_ics["ic"]],
        "x",
        label="hklee_algo_py",
    )
    plt.plot(
        detected_ics_matlab["ic"],
        imu_data["acc_is"].iloc[detected_ics_matlab["ic"]],
        "+",
        label="hklee_algo_matlab",
    )
    plt.xlim(reference_wbs.iloc[3]["start"] - 50, reference_wbs.iloc[3]["end"] + 50)
    plt.legend()
    plt.show()




.. image-sg:: /auto_examples/initial_contacts/images/sphx_glr__03_hklee_algo_001.png
   :alt:  03 hklee algo
   :srcset: /auto_examples/initial_contacts/images/sphx_glr__03_hklee_algo_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 142-146

Evaluation of the algorithm against a reference
--------------------------------------------------
To quantify how the Python output compares to the reference labels, we are providing a range of evaluation functions.
See the :ref:`example on ICD evaluation <icd_evaluation>` for more details.


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.300 seconds)

**Estimated memory usage:**  9 MB


.. _sphx_glr_download_auto_examples_initial_contacts__03_hklee_algo.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: _03_hklee_algo.py <_03_hklee_algo.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: _03_hklee_algo.ipynb <_03_hklee_algo.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_