Improve auto northing

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "res-wind-up"
-version = "0.4.3"
+version = "0.4.4"
 authors = [
     { name = "Alex Clerc", email = "alex.clerc@res-group.com" }
 ]
@@ -69,7 +69,6 @@ examples = [
     'ephem',
     'flaml[automl]',
 ]
-all = ["res-wind-up[dev,examples]"]
 
 [tool.setuptools.packages.find]
 where = ["."]
@@ -127,7 +126,7 @@ module = [
     "geographiclib.geodesic",
     "pandas",
     "pandas.testing",
-    "ruptures",
+    "ruptures.*",
     "scipy.stats",
     "seaborn",
     "utm",

tests/test_math_funcs.py

@@ -4,8 +4,9 @@
 import pandas as pd
 import pytest
 from pandas.testing import assert_series_equal
+from scipy.stats import circmean
 
-from wind_up.circular_math import circ_diff
+from wind_up.circular_math import circ_diff, circ_median
 
 test_circ_diff_data = [
     (0, 0, 0),
@@ -24,6 +25,7 @@
 def test_circ_diff(angle1: float | np.generic, angle2: float | np.generic, expected: float | np.generic) -> None:
     if isinstance(expected, pd.Series):
         assert_series_equal(circ_diff(angle1, angle2), (expected))
+        assert_series_equal(circ_diff(angle1 - 360, angle2 + 360), (expected))
     else:
         assert circ_diff(angle1, angle2) == pytest.approx(expected)
 
@@ -33,3 +35,152 @@ def test_within_bin() -> None:
     d = 242
     dir_bin_width = 10.0
     assert not np.abs(circ_diff(d - 5, d)) < dir_bin_width / 2
+
+
+def circ_median_exact(angles: np.ndarray | list) -> float:
+    """Exact circular median using O(n²) approach for testing.
+
+    In case of a tie, returns the circular mean of the tied angles.
+    """
+
+    angles = np.asarray(angles).flatten()
+    angles = angles[~np.isnan(angles)]
+
+    if len(angles) == 0:
+        return np.nan
+
+    angles = np.mod(angles, 360)
+
+    def _sum_circ_dist(candidate: float) -> float:
+        diffs = circ_diff(angles, candidate)
+        return np.sum(np.abs(diffs))
+
+    distances = np.array([_sum_circ_dist(angle) for angle in angles])
+    min_distance = np.min(distances)
+
+    # Find all angles that have the minimum distance (handle ties)
+    tied_indices = np.where(distances == min_distance)[0]
+
+    if len(tied_indices) == 1:
+        return angles[tied_indices[0]]
+    # Return circular mean of tied angles
+    tied_angles = angles[tied_indices]
+    return circmean(tied_angles, high=360, low=0) % 360
+
+
+test_circ_median_data = [
+    # Simple cases
+    ([0], 0, True),
+    ([0, 20], 10, True),
+    ([0, 15, 20], 15, True),
+    ([350, 0, 15], 0, True),
+    ([170, 181, 190], 181, True),
+    # Edge cases around 0/360 boundary
+    ([355, 0, 15], 0, True),
+    ([350, 351, 9, 10], 0, True),
+    # Symmetric cases
+    ([0, 90, 180, 270], None, True),  # Any answer is valid due to symmetry
+    ([0, 120, 240], None, True),
+    # Single value
+    ([42], 42, True),
+    # Two values
+    ([10, 20], None, True),  # Either 10 or 20 is valid
+    ([350, 10], None, True),  # Could be either
+    # Larger datasets
+    (list(range(10, 351, 1)), 180, True),  # Should be near middle
+    ([i % 360 for i in range(-178, 181, 1)], 1, True),
+    # Test range_360=False
+    ([170, 180, 190], 180, False),
+    ([350, 0, 10], 0, False),
+    ([-10, 0, 10], 0, False),
+]
+
+
+@pytest.mark.parametrize(("angles", "expected", "range_360"), test_circ_median_data)
+def test_circ_median(angles: list, *, expected: float | None, range_360: bool) -> None:
+    result = circ_median(angles, range_360=range_360)
+    if expected is None:
+        assert not np.isnan(result)
+        return
+
+    exact_result = circ_median_exact(angles)
+    exact_result = np.mod(exact_result, 360) if range_360 else np.mod(exact_result + 180, 360) - 180
+
+    # Check that fast and exact methods give similar results
+    # Allow for small differences due to approximation (within 10 degrees)
+    abs_circ_distance = abs(circ_diff(result, exact_result))
+    assert abs_circ_distance < 1e-3, (
+        f"Fast method result {result} differs from exact {exact_result} by {abs_circ_distance} degrees"
+    )
+
+    expected = np.mod(expected, 360) if range_360 else np.mod(expected + 180, 360) - 180
+    abs_circ_distance_expected = abs(circ_diff(result, expected))
+    assert abs_circ_distance_expected < 1e-3, (
+        f"Result {result} differs from expected {expected} by {abs_circ_distance_expected} degrees"
+    )
+
+
+def test_circ_median_with_series() -> None:
+    """Test that pandas Series work correctly."""
+    angles = pd.Series([350, 0, 10, 5])
+    result = circ_median(angles)
+    assert isinstance(result, (float, np.floating))
+    assert result == pytest.approx(2.5)
+
+
+def test_circ_median_with_nan() -> None:
+    """Test that NaN values are handled correctly."""
+    angles = [0, 10, np.nan, 20]
+    result = circ_median(angles)
+    assert not np.isnan(result)
+    assert result == pytest.approx(10)
+
+
+def test_circ_median_all_nan() -> None:
+    """Test that all NaN returns NaN."""
+    angles = [np.nan, np.nan, np.nan]
+    result = circ_median(angles)
+    assert np.isnan(result)
+
+
+@pytest.mark.parametrize("range_360", [True, False])
+def test_circ_median_groupby(*, range_360: bool) -> None:
+    """Test usage with pandas groupby."""
+    df = pd.DataFrame({"group": ["A", "A", "A", "B", "B", "B"], "angle": [350, 0, 15, 170, 181, 195]})
+    result = df.groupby("group")["angle"].apply(lambda x: circ_median(x, range_360=range_360))
+
+    assert len(result) == 2
+    if range_360:
+        assert result["A"] == pytest.approx(0, abs=1e-3)
+        assert result["B"] == pytest.approx(181, abs=1e-3)
+    else:
+        assert result["A"] == pytest.approx(0, abs=1e-3)
+        assert result["B"] == pytest.approx(-179, abs=1e-3)
+
+
+def test_circ_median_performance_comparison() -> None:
+    """Verify that results are consistent between fast and exact methods on larger dataset."""
+    rng = np.random.default_rng(0)
+    # Generate lots of angles near 0 degrees
+    angles = np.concatenate([rng.normal(4, 20, 500) % 360, rng.normal(354, 20, 500) % 360])
+
+    result_fast = circ_median(angles)
+    result_exact = circ_median_exact(angles)
+
+    abs_circ_distance = abs(circ_diff(result_fast, result_exact))
+    assert abs_circ_distance < 1e-1, f"Fast and exact methods differ by {abs_circ_distance} degrees on large dataset"
+
+
+def test_circ_median_range_conversion() -> None:
+    """Test that range_360 parameter works correctly."""
+    angles = [350, 0, 10]
+
+    result_360 = circ_median(angles, range_360=True)
+    assert 0 <= result_360 < 360
+
+    result_180 = circ_median(angles, range_360=False)
+    assert -180 <= result_180 < 180
+
+    # They should represent the same angle
+    abs_circ_distance = abs(circ_diff(result_360, result_180))
+    assert abs_circ_distance < 1e-3

tests/test_northing.py

@@ -3,6 +3,7 @@
 import pandas as pd
 import pytest
 
+from wind_up.circular_math import circ_median
 from wind_up.constants import REANALYSIS_WD_COL
 from wind_up.models import WindUpConfig
 from wind_up.northing import _calc_max_abs_north_errs, apply_northing_corrections
@@ -22,26 +23,26 @@ def test_apply_northing_corrections(test_lsa_t13_config: WindUpConfig) -> None:
     wf_df = _scada_multi_index(wf_df)
     wf_df_after_northing = apply_northing_corrections(wf_df, cfg=cfg, north_ref_wd_col=REANALYSIS_WD_COL, plot_cfg=None)
 
-    median_yaw_before_northing = wf_df.groupby("TurbineName")["YawAngleMean"].median()
-    median_yaw_after_northing = wf_df_after_northing.groupby("TurbineName")["YawAngleMean"].median()
-    assert median_yaw_before_northing.min() == pytest.approx(232.67355346679688)
-    assert median_yaw_before_northing.max() == pytest.approx(232.67355346679688)
-    assert median_yaw_after_northing["LSA_T07"] == pytest.approx(219.69126892089844)
-    assert median_yaw_after_northing["LSA_T09"] == pytest.approx(245.49956512451172)
-    assert median_yaw_after_northing["LSA_T12"] == pytest.approx(177.40547943115234)
-    assert median_yaw_after_northing["LSA_T13"] == pytest.approx(235.22855377197266)
-    assert median_yaw_after_northing["LSA_T14"] == pytest.approx(224.92881774902344)
+    median_yaw_before_northing = wf_df.groupby("TurbineName")["YawAngleMean"].apply(circ_median)
+    median_yaw_after_northing = wf_df_after_northing.groupby("TurbineName")["YawAngleMean"].apply(circ_median)
+    assert median_yaw_before_northing.min() == pytest.approx(237.45018005371094)
+    assert median_yaw_before_northing.max() == pytest.approx(237.45018005371094)
+    assert median_yaw_after_northing["LSA_T07"] == pytest.approx(222.450180)
+    assert median_yaw_after_northing["LSA_T09"] == pytest.approx(253.450180)
+    assert median_yaw_after_northing["LSA_T12"] == pytest.approx(77.712486)
+    assert median_yaw_after_northing["LSA_T13"] == pytest.approx(240.450180)
+    assert median_yaw_after_northing["LSA_T14"] == pytest.approx(228.450180)
 
     abs_north_errs_before_northing = _calc_max_abs_north_errs(
         wf_df, north_ref_wd_col=REANALYSIS_WD_COL, timebase_s=cfg.timebase_s
     )
     abs_north_errs_after_northing = _calc_max_abs_north_errs(
         wf_df_after_northing, north_ref_wd_col=REANALYSIS_WD_COL, timebase_s=cfg.timebase_s
     )
-    assert abs_north_errs_before_northing.min() == pytest.approx(7.88920288085938)
-    assert abs_north_errs_before_northing.max() == pytest.approx(7.88920288085938)
-    assert abs_north_errs_after_northing["LSA_T07"] == pytest.approx(18.400006103515604)
-    assert abs_north_errs_after_northing["LSA_T09"] == pytest.approx(23.889203)
-    assert abs_north_errs_after_northing["LSA_T12"] == pytest.approx(162.918431)
-    assert abs_north_errs_after_northing["LSA_T13"] == pytest.approx(10.889203)
-    assert abs_north_errs_after_northing["LSA_T14"] == pytest.approx(12.400006)
+    assert abs_north_errs_before_northing.min() == pytest.approx(7.911393667045218)
+    assert abs_north_errs_before_northing.max() == pytest.approx(7.911393667045218)
+    assert abs_north_errs_after_northing["LSA_T07"] == pytest.approx(18.402401473162058)
+    assert abs_north_errs_after_northing["LSA_T09"] == pytest.approx(174.10443861846687)
+    assert abs_north_errs_after_northing["LSA_T12"] == pytest.approx(172.59341754958666)
+    assert abs_north_errs_after_northing["LSA_T13"] == pytest.approx(10.894109266368332)
+    assert abs_north_errs_after_northing["LSA_T14"] == pytest.approx(12.41110384502656)

tests/test_optimize_northing.py

@@ -6,6 +6,7 @@
 from pandas.testing import assert_frame_equal
 
 from tests.conftest import TEST_DATA_FLD
+from wind_up.circular_math import circ_median
 from wind_up.constants import RAW_DOWNTIME_S_COL, RAW_POWER_COL, RAW_YAWDIR_COL, TIMESTAMP_COL
 from wind_up.models import WindUpConfig
 from wind_up.optimize_northing import _clip_wtg_north_table, auto_northing_corrections
@@ -105,7 +106,16 @@ def test_clip_wtg_north_table_entry_after_start() -> None:
     assert_frame_equal(actual_wtg_north_table, expected_wtg_north_table)
 
 
-def test_auto_northing_corrections(test_homer_config: WindUpConfig) -> None:
+wind_direction_offsets = [
+    0,
+    343,  # chosen to create lots of 0-360 wraps in the original data
+    290,  # chosen to create lots of 0-360 wraps in the northed data
+]
+
+
+@pytest.mark.slow
+@pytest.mark.parametrize(("wind_direction_offset"), wind_direction_offsets)
+def test_auto_northing_corrections(test_homer_config: WindUpConfig, wind_direction_offset: float) -> None:
     cfg = test_homer_config
     cfg.lt_first_dt_utc_start = pd.Timestamp("2023-07-01 00:00:00", tz="UTC")
     cfg.analysis_last_dt_utc_start = pd.Timestamp("2023-07-31 23:50:00", tz="UTC")
@@ -121,15 +131,25 @@ def test_auto_northing_corrections(test_homer_config: WindUpConfig) -> None:
     wf_df[RAW_YAWDIR_COL] = wf_df["YawAngleMean"]
     wf_df[RAW_DOWNTIME_S_COL] = wf_df["ShutdownDuration"]
 
+    # add wind_direction_offset to direction columns
+    for col in {RAW_YAWDIR_COL, "YawAngleMean", "reanalysis_wd"}:
+        wf_df[col] = (wf_df[col] + wind_direction_offset) % 360
+    if wind_direction_offset != 0:
+        # in this case YawAngleMin and YawAngleMax will be incorrect, so nan them out
+        wf_df["YawAngleMin"] = np.nan
+        wf_df["YawAngleMax"] = np.nan
+
     northed_wf_df = auto_northing_corrections(wf_df, cfg=cfg, plot_cfg=None)
 
-    median_yaw_before_northing = wf_df.groupby("TurbineName", observed=True)["YawAngleMean"].median()
-    median_yaw_after_northing = northed_wf_df.groupby("TurbineName", observed=True)["YawAngleMean"].median()
+    median_yaw_before_northing = wf_df.groupby("TurbineName", observed=True)["YawAngleMean"].apply(circ_median)
+    median_yaw_after_northing = northed_wf_df.groupby("TurbineName", observed=True)["YawAngleMean"].apply(circ_median)
 
-    assert median_yaw_before_northing["HMR_T01"] == pytest.approx(191.0)
-    assert median_yaw_before_northing["HMR_T02"] == pytest.approx(173.0)
-    assert median_yaw_after_northing["HMR_T01"] == pytest.approx(269.6)
-    assert median_yaw_after_northing["HMR_T02"] == pytest.approx(267.6)
+    expected_t1_yaw_after_northing = (290 + wind_direction_offset) % 360
+    expected_t2_yaw_after_northing = (295 + wind_direction_offset) % 360
+    assert median_yaw_before_northing["HMR_T01"] == pytest.approx((343 + wind_direction_offset) % 360)
+    assert median_yaw_before_northing["HMR_T02"] == pytest.approx((173 + wind_direction_offset) % 360)
+    assert median_yaw_after_northing["HMR_T01"] == pytest.approx(expected_t1_yaw_after_northing, abs=1.0)
+    assert median_yaw_after_northing["HMR_T02"] == pytest.approx(expected_t2_yaw_after_northing, abs=1.0)
 
     # try to mess up the yaw angles further and run again
     wf_df[RAW_YAWDIR_COL] = (wf_df[RAW_YAWDIR_COL] + 180) % 360
@@ -149,10 +169,6 @@ def test_auto_northing_corrections(test_homer_config: WindUpConfig) -> None:
 
     northed_wf_df = auto_northing_corrections(wf_df, cfg=cfg, plot_cfg=None)
 
-    median_yaw_before_northing = wf_df.groupby("TurbineName", observed=True)["YawAngleMean"].median()
-    median_yaw_after_northing = northed_wf_df.groupby("TurbineName", observed=True)["YawAngleMean"].median()
-
-    assert median_yaw_before_northing["HMR_T01"] == pytest.approx(178.0)
-    assert median_yaw_before_northing["HMR_T02"] == pytest.approx(206.0)
-    assert median_yaw_after_northing["HMR_T01"] == pytest.approx(269.2)
-    assert median_yaw_after_northing["HMR_T02"] == pytest.approx(269.2)
+    median_yaw_after_northing = northed_wf_df.groupby("TurbineName", observed=True)["YawAngleMean"].apply(circ_median)
+    assert median_yaw_after_northing["HMR_T01"] == pytest.approx(expected_t1_yaw_after_northing, abs=1.5)
+    assert median_yaw_after_northing["HMR_T02"] == pytest.approx(expected_t2_yaw_after_northing, abs=1.5)