fourier series: more consistent with standard definition (MESMER-grou…

…p#512)

* fourier series: more consistent with standard definition

* changelog

CHANGELOG.rst

@@ -147,10 +147,11 @@ Harmonic model
 - de-duplicate the expression of months in their harmonic form (`#415 <https://github.com/MESMER-group/mesmer/pull/415>`_)
   move creation of the month array to the deepest level (`#487 <https://github.com/MESMER-group/mesmer/pull/487>`_).
 - fix indexing of harmonic model coefficients (`#415 <https://github.com/MESMER-group/mesmer/pull/415>`_)
--  Refactor variable names, small code improvements, fixes and clean docstring (
-   `#415 <https://github.com/MESMER-group/mesmer/pull/415>`_,
-   `#424 <https://github.com/MESMER-group/mesmer/pull/424>`_, and
-   `#433 <https://github.com/MESMER-group/mesmer/pull/433>`_)
+-  Refactor variable names, small code improvements, fixes and clean docstring
+   (`#415 <https://github.com/MESMER-group/mesmer/pull/415>`_,
+   `#424 <https://github.com/MESMER-group/mesmer/pull/424>`_,
+   `#433 <https://github.com/MESMER-group/mesmer/pull/433>`_, and
+   `#512 <https://github.com/MESMER-group/mesmer/pull/512>`_)
 - add tests (
   `#431 <https://github.com/MESMER-group/mesmer/pull/431>`_, and
   `#458 <https://github.com/MESMER-group/mesmer/pull/458>`_)

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,

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,
         ]
     )