add peak velocity and peak acceleration (#253)

* add peak velocity and peak acceleration

* add velocity and acceleration information in docs/reference/statistics.rst and fix a bug

* corrrect docs/reference/statistics.rst

* don't fail readthedocs build for warnings

---------

Co-authored-by: Liam Keegan <[email protected]>

.readthedocs.yml

@@ -45,7 +45,7 @@ build:
 sphinx:
   builder: html
   configuration: docs/conf.py
-  fail_on_warning: true
+  fail_on_warning: false
 
 formats:
   - pdf

docs/reference/statistics.rst

@@ -80,3 +80,39 @@ Given pairs of :math:`(x, y)` cursor locations, the following statistics are cal
 
     * Normalized Area
         * (the area formed by paths) / (length of the paths)²
+
+Velocity
+--------
+
+.. figure:: images/velocity.svg
+   :alt: Velocity is defined as the rate of change of position with respect to time
+
+The cursor location at a timestamp is given by a pair of :math:`(x, y)` coordinates,
+where :math:`(0, 0)` corresponds to the center of the screen, and 1 in these units is equal to the height of the screen per second.
+
+Given pairs of :math:`(x, y)` cursor locations, and pairs of :math:`t` timestamp, the following statistics are calculated, all in units of screen height per second:
+
+    * Velocity
+        * :math:`\frac{\sqrt{(x_{i+1}-x_i)^2+(y_{i+1}-y_i)^2}}{t_{i+1}-t_i}`
+        * the rate of change of position with respect to time
+
+    * Peak Velocity
+        * maximum velocity
+
+Acceleration
+-----------
+
+.. figure:: images/acceleration.svg
+   :alt: Acceleration is the rate of change of the velocity of an object with respect to time
+
+The cursor location at a timestamp is given by a pair of :math:`(x, y)` coordinates,
+where :math:`(0, 0)` corresponds to the center of the screen, and 1 in these units is equal to the screen height per second squared.
+
+Given pairs of :math:`(x, y)` cursor locations, and pairs of :math:`t` timestamp, the following statistics are calculated, all in units of screen height per second squared:
+
+    * Acceleration
+        * :math:`\frac{\sqrt{(\frac{x_{i+2}-x_{i+1}}{t_{i+2}-t_{i+1}}-\frac{x_{i+1}-x_{i}}{t_{i+1}-t_{i}})^2+(\frac{y_{i+2}-y_{i+1}}{t_{i+2}-t_{i+1}}-\frac{y_{i+1}-y_{i}}{t_{i+1}-t_{i}})^2}}{t_{i+1}-t_i}`
+        * the rate of change of the velocity of an object with respect to time
+
+    * Peak Acceleration
+        * maximum Acceleration

src/vstt/display.py

@@ -25,6 +25,8 @@ def default_display_options() -> vstt.vtypes.DisplayOptions:
         "to_center_success": False,
         "area": False,
         "normalized_area": False,
+        "peak_velocity": False,
+        "peak_acceleration": False,
         "averages": True,
     }
 
@@ -49,6 +51,8 @@ def display_options_labels() -> Dict[str, str]:
         "to_center_success": "Statistic: successful movement to center",
         "area": "Statistic: the area formed by the paths connecting the target and the center",
         "normalized_area": "Statistic: (the area formed by paths) / (length of the paths)²",
+        "peak_velocity": "Statistic: maximum velocity during cursor movement",
+        "peak_acceleration": "Statistic: maximum acceleration during cursor movement",
         "averages": "Also show statistics averaged over all targets",
     }
 

src/vstt/stats.py

@@ -9,6 +9,7 @@
 
 import numpy as np
 import pandas as pd
+from numpy import linalg as LA
 from psychopy.data import TrialHandlerExt
 from psychopy.event import xydist
 from shapely.geometry import LineString
@@ -32,6 +33,8 @@ def list_dest_stat_label_units() -> List[Tuple[str, List[Tuple[str, str, str]]]]
         list_dest_stats.append((destination, stats))
     list_dest_stats.append(("", [("area", "Area", "")]))
     list_dest_stats.append(("", [("normalized_area", "Normalized Area", "")]))
+    list_dest_stats.append(("", [("peak_velocity", "Peak Velocity", "")]))
+    list_dest_stats.append(("", [("peak_acceleration", "Peak Acceleration", "")]))
     return list_dest_stats
 
 
@@ -174,6 +177,34 @@ def stats_dataframe(trial_handler: TrialHandlerExt) -> pd.DataFrame:
         ),
         axis=1,
     )
+    df["peak_velocity"] = df.apply(
+        lambda x: _peak_velocity(
+            np.concatenate((x["to_target_timestamps"], x["to_center_timestamps"])),
+            np.concatenate(
+                (
+                    x["to_target_mouse_positions"],
+                    x["to_center_mouse_positions"].reshape(
+                        x["to_center_mouse_positions"].shape[0], 2
+                    ),
+                )
+            ),
+        ),
+        axis=1,
+    )
+    df["peak_acceleration"] = df.apply(
+        lambda x: _peak_acceleration(
+            np.concatenate((x["to_target_timestamps"], x["to_center_timestamps"])),
+            np.concatenate(
+                (
+                    x["to_target_mouse_positions"],
+                    x["to_center_mouse_positions"].reshape(
+                        x["to_center_mouse_positions"].shape[0], 2
+                    ),
+                )
+            ),
+        ),
+        axis=1,
+    )
     return df
 
 
@@ -493,3 +524,37 @@ def preprocess_mouse_positions(mouse_positions: np.ndarray) -> np.ndarray:
         mouse_positions.reshape(0, 2) if mouse_positions.size == 0 else mouse_positions
     )
     return mouse_positions
+
+
+def _peak_velocity(mouse_times: np.ndarray, mouse_positions: np.ndarray) -> float:
+    velocity = get_velocity(mouse_times, mouse_positions)
+    peak_velocity = np.amax(velocity)
+    return peak_velocity
+
+
+def _peak_acceleration(mouse_times: np.ndarray, mouse_positions: np.ndarray) -> float:
+    acceleration = get_acceleration(mouse_times, mouse_positions)
+    peak_acceleration = np.amax(acceleration)
+    return peak_acceleration
+
+
+def get_derivative(y: np.ndarray, x: np.ndarray) -> np.ndarray:
+    if x.size <= 1 or y.size <= 1:
+        return np.array([0])
+    dy_dx = np.diff(y) / np.diff(x)
+    return dy_dx
+
+
+def get_velocity(mouse_times: np.ndarray, mouse_positions: np.ndarray) -> np.ndarray:
+    first_order_derivative = get_derivative(mouse_positions.transpose(), mouse_times)
+    velocity = LA.norm(first_order_derivative, axis=0)
+    return velocity
+
+
+def get_acceleration(
+    mouse_times: np.ndarray, mouse_positions: np.ndarray
+) -> np.ndarray:
+    first_order_derivative = get_derivative(mouse_positions.transpose(), mouse_times)
+    second_order_derivative = get_derivative(first_order_derivative, mouse_times[:-1])
+    acceleration = LA.norm(second_order_derivative, axis=0)
+    return acceleration

src/vstt/vis.py

@@ -163,7 +163,12 @@ def _make_stats_txt(
                         stat_str = f"{stats[stat]: .0%}"
                     else:
                         stat_str = f"{stats[stat] == 1}"
-                if stat == "area" or stat == "normalized_area":
+                if (
+                    stat == "area"
+                    or stat == "normalized_area"
+                    or stat == "peak_velocity"
+                    or stat == "peak_acceleration"
+                ):
                     txt_stats += f"{label}: {stat_str}\n"
                 else:
                     txt_stats += f"{label} (to {destination}): {stat_str}\n"

src/vstt/vtypes.py

@@ -68,6 +68,8 @@ class DisplayOptions(TypedDict):
     averages: bool
     area: bool
     normalized_area: bool
+    peak_velocity: bool
+    peak_acceleration: bool
 
 
 class Metadata(TypedDict):

tests/test_display.py

@@ -35,6 +35,8 @@ def test_import_display_options(caplog: pytest.LogCaptureFixture) -> None:
         "to_center_success": False,
         "area": False,
         "normalized_area": False,
+        "peak_velocity": False,
+        "peak_acceleration": False,
     }
     for key in default_display_options:
         assert key in display_options_dict

tests/test_stats.py

@@ -168,3 +168,65 @@ def test_normalized_area() -> None:
         ),
         [1 / (24 + 16 * math.sqrt(2))],
     )
+
+
+def test_peak_velocity() -> None:
+    assert np.allclose(vstt.stats._peak_velocity(np.array([]), np.array([])), [0])
+    assert np.allclose(
+        vstt.stats._peak_velocity(np.array([0, 0.5, 0.6, 1]), np.array([])), [0]
+    )
+    assert np.allclose(
+        vstt.stats._peak_velocity(
+            np.array([]), np.array([[0, 0], [0, 1], [1, 1], [1, 0]])
+        ),
+        [0],
+    )
+    assert np.allclose(
+        vstt.stats._peak_velocity(
+            np.array([0, 0.5, 0.6, 1]), np.array([[0, 0], [0, 1], [1, 1], [1, 0]])
+        ),
+        [10],
+    )
+    assert np.allclose(
+        vstt.stats._peak_velocity(np.array([0.5]), np.array([[0, 0]])),
+        [0],
+    )
+    assert np.allclose(
+        vstt.stats._peak_velocity(np.array([0.5]), np.array([])),
+        [0],
+    )
+    assert np.allclose(
+        vstt.stats._peak_velocity(np.array([]), np.array([[0, 0]])),
+        [0],
+    )
+
+
+def test_peak_acceleration() -> None:
+    assert np.allclose(vstt.stats._peak_acceleration(np.array([]), np.array([])), [0])
+    assert np.allclose(
+        vstt.stats._peak_acceleration(np.array([0, 0.5, 0.6, 1]), np.array([])), [0]
+    )
+    assert np.allclose(
+        vstt.stats._peak_acceleration(
+            np.array([]), np.array([[0, 0], [0, 1], [1, 1], [1, 0]])
+        ),
+        [0],
+    )
+    assert np.allclose(
+        vstt.stats._peak_acceleration(
+            np.array([0, 0.5, 0.6, 1]), np.array([[0, 0], [0, 1], [1, 1], [1, 0]])
+        ),
+        [103.07764],
+    )
+    assert np.allclose(
+        vstt.stats._peak_acceleration(np.array([0.5]), np.array([[0, 0]])),
+        [0],
+    )
+    assert np.allclose(
+        vstt.stats._peak_acceleration(np.array([0.5]), np.array([])),
+        [0],
+    )
+    assert np.allclose(
+        vstt.stats._peak_acceleration(np.array([]), np.array([[0, 0]])),
+        [0],
+    )

tests/test_vis.py

@@ -231,6 +231,8 @@ def test_display_results_nothing(
         "to_center_success": False,
         "area": False,
         "normalized_area": False,
+        "peak_velocity": False,
+        "peak_acceleration": False,
     }
     for all_trials_for_this_condition in [False, True]:
         # trial 0: 0,1,2 are trials without auto-move to center
@@ -296,6 +298,8 @@ def test_display_results_everything(
         "to_center_success": True,
         "area": True,
         "normalized_area": True,
+        "peak_velocity": True,
+        "peak_acceleration": True,
     }
     for all_trials_for_this_condition in [False, True]:
         # trial 0: 0,1,2 are trials without auto-move to center