Skip to content

Commit

Permalink
add: utils/rotate_data()
Browse files Browse the repository at this point in the history
  • Loading branch information
huangziwei committed Jan 29, 2025
1 parent d5cc0a3 commit 39d9f9e
Show file tree
Hide file tree
Showing 5 changed files with 147 additions and 57 deletions.
17 changes: 9 additions & 8 deletions pycircstat2/base.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from typing import Union
from typing import Optional, Union

import numpy as np

Expand All @@ -15,7 +15,7 @@
circ_std,
)
from .hypothesis import rayleigh_test
from .utils import data2rad, rad2data, significance_code
from .utils import data2rad, rad2data, rotate_data, significance_code
from .visualization import circ_plot

__names__ = ["Circular", "Axial"]
Expand Down Expand Up @@ -146,13 +146,14 @@ class Circular:
def __init__(
self,
data: Union[np.ndarray, list], # angle
w: Union[np.ndarray, list, None] = None, # frequency
bins: Union[int, np.ndarray, None] = None,
w: Optional[Union[np.ndarray, list]] = None, # frequency
bins: Optional[Union[int, np.ndarray]] = None,
unit: str = "degree",
n_intervals: Union[
int, float, None
] = None, # number of intervals in the full cycle
n_intervals: Optional[Union[
int, float
]] = None, # number of intervals in the full cycle
n_clusters_max: int = 1, # number of clusters to be tested for mixture of von Mises
rotate: Optional[float] = None, # in rad
**kwargs,
):
# meta
Expand Down Expand Up @@ -184,7 +185,7 @@ def __init__(

# data
self.data = data = np.array(data) if isinstance(data, list) else data
self.alpha = alpha = data2rad(data, n_intervals)
self.alpha = alpha = data2rad(data, n_intervals) if rotate is None else rotate_data(data2rad(data, n_intervals), rotate, unit="radian")

# data preprocessing
if bins is None:
Expand Down
96 changes: 52 additions & 44 deletions pycircstat2/descriptive.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from typing import Tuple, Union
from typing import Optional, Tuple, Union

import numpy as np
from scipy.stats import chi2, norm
Expand All @@ -7,10 +7,10 @@


def circ_r(
alpha: Union[np.ndarray, None] = None,
w: Union[np.ndarray, None] = None,
Cbar: Union[float, None] = None,
Sbar: Union[float, None] = None,
alpha: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
Cbar: Optional[float] = None,
Sbar: Optional[float] = None,
) -> float:
r"""
Circular mean resultant vector length (r).
Expand Down Expand Up @@ -54,7 +54,7 @@ def circ_r(

def circ_mean(
alpha: np.ndarray,
w: Union[np.ndarray, None] = None,
w: Optional[np.ndarray] = None,
) -> Union[np.ndarray, float]:
r"""
Circular mean (m).
Expand Down Expand Up @@ -106,7 +106,7 @@ def circ_mean(

def circ_mean_and_r(
alpha: np.ndarray,
w: Union[np.ndarray, None] = None,
w: Optional[np.ndarray] = None,
) -> Tuple[Union[float, np.ndarray], float]:
"""
Circular mean (m) and resultant vector length (r).
Expand Down Expand Up @@ -148,8 +148,8 @@ def circ_mean_and_r(

def circ_mean_and_r_of_means(
circs: Union[list, None] = None,
ms: Union[np.ndarray, None] = None,
rs: Union[np.ndarray, None] = None,
ms: Optional[np.ndarray] = None,
rs: Optional[np.ndarray] = None,
) -> Tuple[float, float]:
"""The Mean of a set of Mean Angles
Expand Down Expand Up @@ -194,7 +194,7 @@ def circ_mean_and_r_of_means(

def circ_moment(
alpha: np.ndarray,
w: Union[np.ndarray, None] = None,
w: Optional[np.ndarray] = None,
p: int = 1,
mean: Union[float, np.ndarray, None] = None,
centered: bool = False,
Expand Down Expand Up @@ -241,7 +241,7 @@ def circ_moment(

def circ_dispersion(
alpha: np.ndarray,
w: Union[np.ndarray, None] = None,
w: Optional[np.ndarray] = None,
mean=None,
) -> float:
r"""
Expand Down Expand Up @@ -281,7 +281,7 @@ def circ_dispersion(
return dispersion


def circ_skewness(alpha: np.ndarray, w: Union[np.ndarray, None] = None) -> float:
def circ_skewness(alpha: np.ndarray, w: Optional[np.ndarray] = None) -> float:
r"""
Circular skewness, as defined by Equation 2.29 (Fisher, 1993):
Expand Down Expand Up @@ -317,7 +317,7 @@ def circ_skewness(alpha: np.ndarray, w: Union[np.ndarray, None] = None) -> float
return skewness


def circ_kurtosis(alpha: np.ndarray, w: Union[np.ndarray, None] = None) -> float:
def circ_kurtosis(alpha: np.ndarray, w: Optional[np.ndarray] = None) -> float:
r"""
Circular kurtosis, as defined by Equation 2.30 (Fisher, 1993):
Expand Down Expand Up @@ -354,10 +354,10 @@ def circ_kurtosis(alpha: np.ndarray, w: Union[np.ndarray, None] = None) -> float


def angular_var(
alpha: Union[np.ndarray, None] = None,
w: Union[np.ndarray, None] = None,
r: Union[float, None] = None,
bin_size: Union[float, None] = None,
alpha: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
r: Optional[float] = None,
bin_size: Optional[float] = None,
) -> float:
r"""
Angular variance
Expand Down Expand Up @@ -391,10 +391,10 @@ def angular_var(


def angular_std(
alpha: Union[np.ndarray, None] = None,
w: Union[np.ndarray, None] = None,
r: Union[float, None] = None,
bin_size: Union[float, None] = None,
alpha: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
r: Optional[float] = None,
bin_size: Optional[float] = None,
) -> float:
r"""
Angular (standard) deviation
Expand Down Expand Up @@ -430,10 +430,10 @@ def angular_std(


def circ_var(
alpha: Union[np.ndarray, None] = None,
w: Union[np.ndarray, None] = None,
r: Union[float, None] = None,
bin_size: Union[float, None] = None,
alpha: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
r: Optional[float] = None,
bin_size: Optional[float] = None,
) -> float:
r"""
Circular variance
Expand Down Expand Up @@ -494,10 +494,10 @@ def circ_var(


def circ_std(
alpha: Union[np.ndarray, None] = None,
w: Union[np.ndarray, None] = None,
r: Union[float, None] = None,
bin_size: Union[float, None] = None,
alpha: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
r: Optional[float] = None,
bin_size: Optional[float] = None,
) -> tuple:
r"""
Circular standard deviation (s).
Expand Down Expand Up @@ -535,7 +535,7 @@ def circ_std(

def circ_median(
alpha: np.ndarray,
w: Union[np.ndarray, None] = None,
w: Optional[np.ndarray] = None,
method: str = "deviation",
return_average: bool = True,
average_method: str = "all",
Expand Down Expand Up @@ -817,10 +817,10 @@ def circ_mean_deviation(


def circ_mean_ci(
alpha: Union[np.ndarray, None] = None,
w: Union[np.ndarray, None] = None,
mean: Union[float, None] = None,
r: Union[float, None] = None,
alpha: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
mean: Optional[float] = None,
r: Optional[float] = None,
n: Union[int, None] = None,
ci: float = 0.95,
method: str = "approximate",
Expand Down Expand Up @@ -921,8 +921,8 @@ def circ_mean_ci(

def _circ_mean_ci_dispersion(
alpha: np.ndarray,
w: Union[np.ndarray, None] = None,
mean: Union[float, None] = None,
w: Optional[np.ndarray] = None,
mean: Optional[float] = None,
ci: float = 0.95,
) -> tuple[float, float]:
r"""Confidence intervals based on circular dispersion.
Expand Down Expand Up @@ -978,10 +978,10 @@ def _circ_mean_ci_dispersion(


def _circ_mean_ci_approximate(
alpha: Union[np.ndarray, None] = None,
w: Union[np.ndarray, None] = None,
mean: Union[float, None] = None,
r: Union[float, None] = None,
alpha: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
mean: Optional[float] = None,
r: Optional[float] = None,
n: Union[int, None] = None,
ci: float = 0.95,
) -> tuple:
Expand Down Expand Up @@ -1146,8 +1146,8 @@ def _circ_mean_resample(alpha, z0, v0):

def circ_median_ci(
median: float = None,
alpha: Union[np.ndarray, None] = None,
w: Union[np.ndarray, None] = None,
alpha: Optional[np.ndarray] = None,
w: Optional[np.ndarray] = None,
method: str = "deviation",
ci: float = 0.95,
) -> tuple:
Expand Down Expand Up @@ -1201,10 +1201,18 @@ def circ_median_ci(

offset = int(1 + np.floor(0.5 * np.sqrt(n) * z)) # fisher:eq(4.19)

idx_median = np.where(alpha.round(5) < median.round(5))[0][-1]
# idx_median = np.where(alpha.round(5) < np.round(median, 5))[0][-1]
arr = np.where(alpha.round(5) < np.round(median, 5))[0]
if len(arr) == 0:
# That means median is smaller than alpha[0] (to 5 decimals).
# In a circular sense, the “closest index below” is alpha[-1].
idx_median = len(alpha) - 1
else:
idx_median = arr[-1]

idx_lb = idx_median - offset + 1
idx_ub = idx_median + offset
if median.round(5) in alpha.round(5): # don't count the median per se
if np.round(median, 5) in alpha.round(5): # don't count the median per se
idx_ub += 1

if idx_ub > n:
Expand Down
41 changes: 40 additions & 1 deletion pycircstat2/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,7 @@ def load_data(
return csv_data


def is_within_circular_range(value, lb, ub):
def is_within_circular_range(value: float, lb: float, ub: float) -> bool:
"""
Check if a value lies within the circular range [lb, ub].
Expand Down Expand Up @@ -188,3 +188,42 @@ def is_within_circular_range(value, lb, ub):
else:
# Wrapping range
return value >= lb or value <= ub

def rotate_data(alpha: np.ndarray, angle: float, unit: str = "radian") -> np.ndarray:
"""
Rotate a circular dataset by a given angle, supporting degrees, radians, and hours.
Parameters
----------
alpha : np.ndarray
Angles in the specified unit.
angle : float
Rotation angle in the specified unit.
unit : str, optional
Unit of measurement ("degree", "radian", or "hour"). Default is "radian".
Returns
-------
rotated_alpha : np.ndarray
Rotated angles, normalized within the unit's full cycle.
"""
if unit == "degree":
n_intervals = 360
elif unit == "radian":
n_intervals = 2 * np.pi
elif unit == "hour":
n_intervals = 24
else:
raise ValueError("Unit must be 'degree', 'radian', or 'hour'.")

# Convert to radians for consistent computation
alpha_rad = data2rad(alpha, k=n_intervals)
angle_rad = data2rad(angle, k=n_intervals)

# Perform rotation and normalize in radians
rotated_alpha_rad = angmod(alpha_rad + angle_rad, bounds=[0, 2 * np.pi])

# Convert back to the original unit
rotated_alpha = rad2data(rotated_alpha_rad, k=n_intervals)

return rotated_alpha
4 changes: 1 addition & 3 deletions pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,7 @@ readme = {file = "README.md", content-type = "text/markdown"}

[project.optional-dependencies]
dev = [
"black",
"isort",
"jupyter",
"ruff",
"pytest",
]

Expand Down
46 changes: 45 additions & 1 deletion tests/test_utils.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,13 @@
import numpy as np

from pycircstat2.utils import angmod, angular_distance, data2rad, rad2data, time2float
from pycircstat2.utils import (
angmod,
angular_distance,
data2rad,
rad2data,
rotate_data,
time2float,
)


def test_data2rad():
Expand Down Expand Up @@ -84,3 +91,40 @@ def test_angular_distance():
np.rad2deg(angular_distance(a=np.deg2rad(190), b=np.deg2rad(5))).round(3),
175,
)

def test_rotate_data():
# Test case 1: Rotating angles in radians
alpha = np.array([0, np.pi / 2, np.pi, 3 * np.pi / 2])
rotated = rotate_data(alpha, np.pi / 4, unit="radian")
expected = np.array([np.pi / 4, 3 * np.pi / 4, 5 * np.pi / 4, 7 * np.pi / 4])
np.testing.assert_allclose(rotated, expected, atol=1e-6)

# Test case 2: Rotating angles in degrees
alpha_deg = np.array([0., 90., 180., 270.])
rotated_deg = rotate_data(alpha_deg, 45, unit="degree")
expected_deg = np.array([45., 135., 225., 315.])
np.testing.assert_allclose(rotated_deg, expected_deg, atol=1e-6)

# Test case 3: Rotating angles in hours
alpha_hour = np.array([0., 6., 12., 18.])
rotated_hour = rotate_data(alpha_hour, 3, unit="hour")
expected_hour = np.array([3., 9., 15., 21.])
np.testing.assert_allclose(rotated_hour, expected_hour, atol=1e-6)

# Test case 4: Rotation that wraps around the cycle
rotated_wrap = rotate_data(np.array([350.]), 20, unit="degree")
expected_wrap = np.array([10])
np.testing.assert_allclose(rotated_wrap, expected_wrap, atol=1e-6)

rotated_wrap_hour = rotate_data(np.array([22.]), 5, unit="hour")
expected_wrap_hour = np.array([3.])
np.testing.assert_allclose(rotated_wrap_hour, expected_wrap_hour, atol=1e-6)

# Test case 5: Rotating an empty array should return an empty array
np.testing.assert_array_equal(rotate_data(np.array([]), 45., unit="degree"), np.array([]))

# Test case 6: Rotating by 0 should return the same values
alpha_no_rotation = np.array([10., 50., 120.])
np.testing.assert_allclose(
rotate_data(alpha_no_rotation, 0, unit="degree"), alpha_no_rotation
)

0 comments on commit 39d9f9e

Please sign in to comment.