#36: Fix incorrect wave identifier mapping in harmonic analysis for s…

…parse tables of constituents.

include/fes/wave/table.hpp

@@ -143,13 +143,22 @@ class Table {
   /// 2N2, N2, M2, K2.
   auto admittance() -> void;
 
-  /// Get the wave properties
+  /// Get the wave properties from its index
+  inline auto operator[](const size_t index) const
+      -> const std::shared_ptr<Wave>& {
+    if (index >= wave_identifiers_.size()) {
+      throw std::out_of_range("index out of range");
+    }
+    return waves_[wave_identifiers_[index]];
+  }
+
+  /// Get the wave properties from its identifier
   constexpr auto operator[](const Constituent ident) const
       -> const std::shared_ptr<Wave>& {
     return ((*this).*getter_)(ident);
   }
 
-  /// Get the wave properties
+  /// Get the wave properties from its identifier
   constexpr auto operator[](const Constituent ident) -> std::shared_ptr<Wave>& {
     return std::remove_const_t<std::shared_ptr<Wave>&>(
         ((*this).*getter_)(ident));
@@ -323,6 +332,11 @@ class Table {
   /// Wave table
   std::vector<std::shared_ptr<Wave>> waves_{};
 
+  /// An array that maps linear indices (0, 1, 2, 3, ...) to the wave
+  /// identifiers defined in the table. If the table is complete, this mapping
+  /// is an identity mapping {0:0, 1:1, 2:2, ...}.
+  std::vector<size_t> wave_identifiers_{};
+
   /// Get a wave from the table
   ///
   /// @param[in] ident Wave identifier

src/library/wave/table.cpp

@@ -205,6 +205,14 @@ Table::Table(const std::vector<std::string>& waves) {
                 : create_sparse_table(known_constituents, waves, waves_);
   getter_ =
       size() == waves_.size() ? &Table::direct_access : &Table::sparse_access;
+
+  // Fill the index between to have a direct access to the wave
+  wave_identifiers_.reserve(waves_.size());
+  for (const auto& item : waves_) {
+    if (item) {
+      wave_identifiers_.emplace_back(static_cast<size_t>(item->ident()));
+    }
+  }
 }
 
 Table::Table(const std::vector<Constituent>& waves) {
@@ -386,7 +394,7 @@ auto Table::tide_from_tide_series(
     wt.compute_nodal_corrections(angles);
 
     for (size_t jx = 0; jx < wt.size(); ++jx) {
-      const auto& item = wt[static_cast<Constituent>(jx)];
+      const auto& item = wt[jx];
       double phi = item->vu();
 
       tide += item->f() * (wave(jx).real() * std::cos(phi) +
@@ -417,7 +425,7 @@ auto Table::tide_from_mapping(const double epoch, const uint16_t leap_seconds,
 
     for (auto ix = start; ix < end; ++ix) {
       for (size_t jx = 0; jx < wt.size(); ++jx) {
-        const auto& item = wt[static_cast<Constituent>(jx)];
+        const auto& item = wt[jx];
         double phi = item->vu();
 
         result(ix, jx) += item->f() * (wave(jx, ix).real() * std::cos(phi) +
@@ -450,8 +458,8 @@ auto Table::compute_nodal_modulations(
     angles.update(epoch(ix), leap_seconds(ix));
     wt.compute_nodal_corrections(angles);
 
-    for (std::size_t jx = 0; jx < wt.size(); ++jx) {
-      const auto& wave = wt[static_cast<Constituent>(jx)];
+    for (size_t jx = 0; jx < wt.size(); ++jx) {
+      const auto& wave = wt[jx];
       f(jx, ix) = wave->f();
       vu(jx, ix) = wave->vu();
     }

tests/python/core/test_wave_table.py

@@ -6,6 +6,7 @@
 import pathlib
 
 import netCDF4
+import numpy
 from pyfes import core
 from pyfes.leap_seconds import get_leap_seconds
 import pytest
@@ -73,3 +74,24 @@ def test_harmonic_analysis():
 
     assert delta.mean(), pytest.approx(0, rel=1e-16)
     assert delta.std(), pytest.approx(0, rel=1e-12)
+
+
+def test_harmonic_analysis_with_empty_table():
+    time = numpy.arange(numpy.datetime64('2018-01-01'),
+                        numpy.datetime64('2018-03-01'),
+                        3600,
+                        dtype='datetime64[s]')
+    h = numpy.random.rand(time.shape[0])
+    leap_seconds = get_leap_seconds(time)
+
+    wt = core.WaveTable(['M2', 'S2', 'N2', 'K1', 'O1', 'Q1'])
+    w = wt.harmonic_analysis(h,
+                             *wt.compute_nodal_modulations(time, leap_seconds))
+    assert numpy.all(
+        ~numpy.isnan(wt.tide_from_tide_series(time, leap_seconds, w)))
+
+    wt = core.WaveTable()
+    w = wt.harmonic_analysis(h,
+                             *wt.compute_nodal_modulations(time, leap_seconds))
+    assert numpy.all(
+        ~numpy.isnan(wt.tide_from_tide_series(time, leap_seconds, w)))