From 444c682b23b893a876adf9bf651b40c4c108b7d2 Mon Sep 17 00:00:00 2001 From: Mathias Hauser Date: Mon, 9 Sep 2024 11:09:16 +0200 Subject: [PATCH] fourier series: more consistent with standard definition (#512) * fourier series: more consistent with standard definition * changelog --- CHANGELOG.rst | 9 ++++---- mesmer/stats/_harmonic_model.py | 7 ++++--- tests/unit/test_harmonic_model.py | 34 +++++++++++++++---------------- 3 files changed, 26 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 54bde5ef..960b67c8 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -147,10 +147,11 @@ Harmonic model - de-duplicate the expression of months in their harmonic form (`#415 `_) move creation of the month array to the deepest level (`#487 `_). - fix indexing of harmonic model coefficients (`#415 `_) -- Refactor variable names, small code improvements, fixes and clean docstring ( - `#415 `_, - `#424 `_, and - `#433 `_) +- Refactor variable names, small code improvements, fixes and clean docstring + (`#415 `_, + `#424 `_, + `#433 `_, and + `#512 `_) - add tests ( `#431 `_, and `#458 `_) diff --git a/mesmer/stats/_harmonic_model.py b/mesmer/stats/_harmonic_model.py index 0b7a8070..8358b724 100644 --- a/mesmer/stats/_harmonic_model.py +++ b/mesmer/stats/_harmonic_model.py @@ -37,14 +37,15 @@ def _generate_fourier_series_np(yearly_predictor, coeffs): """ order = int(coeffs.size / 4) n_years = yearly_predictor.size // 12 - months = np.tile(np.arange(1, 13), n_years) + # NOTE: months from 0..11 for consistency with standard fourier series and fft + months = np.tile(np.arange(12), n_years) seasonal_cycle = np.nansum( [ (coeffs[idx * 4] * yearly_predictor + coeffs[idx * 4 + 1]) - * np.sin(np.pi * i * (months) / 6) + * np.cos(2 * np.pi * i * months / 12) + (coeffs[idx * 4 + 2] * yearly_predictor + coeffs[idx * 4 + 3]) - * np.cos(np.pi * i * (months) / 6) + * np.sin(2 * np.pi * i * months / 12) for idx, i in enumerate(range(1, order + 1)) ], axis=0, diff --git a/tests/unit/test_harmonic_model.py b/tests/unit/test_harmonic_model.py index bdf0850f..bf457c10 100644 --- a/tests/unit/test_harmonic_model.py +++ b/tests/unit/test_harmonic_model.py @@ -20,21 +20,21 @@ def test_generate_fourier_series_np(): n_months = n_years * 12 yearly_predictor = np.ones(n_months) - months = np.tile(np.arange(1, 13), n_years) + months = np.tile(np.arange(12), n_years) + alpha = 2 * np.pi * months / 12 - coeffs = np.array([0, -1, 0, -2]) + coeffs = np.array([0, -2, 0, -1]) # dummy yearly cycle - expected = coeffs[1] * np.sin(2 * np.pi * (months) / 12) + coeffs[3] * np.cos( - 2 * np.pi * (months) / 12 - ) + expected = coeffs[1] * np.cos(alpha) + coeffs[3] * np.sin(alpha) result = _generate_fourier_series_np(yearly_predictor, coeffs) np.testing.assert_equal(result, expected) - result = _generate_fourier_series_np(yearly_predictor, np.array([3.14, -1, 1, -2])) - expected += 3.14 * np.sin(np.pi * months / 6) + 1 * np.cos(np.pi * months / 6) + coeffs = np.array([1, -2, 3.14, -1]) + result = _generate_fourier_series_np(yearly_predictor, coeffs) + expected += 1 * np.cos(alpha) + 3.14 * np.sin(alpha) np.testing.assert_allclose(result, expected, atol=1e-10) @@ -171,18 +171,18 @@ def test_fit_harmonic_model(): # compare numerically one cell of one year expected = np.array( [ - 9.970548, + 7.324277, 9.966644, - 7.325875, - 2.755833, - -2.518943, - -7.085081, - -9.719088, + 9.972146, + 7.33931, + 2.7736, + -2.501604, + -7.072816, -9.715184, - -7.074415, - -2.504373, - 2.770403, - 7.336541, + -9.720686, + -7.087849, + -2.52214, + 2.753065, ] )