Fixed bug in building convolution windows

- Convolution windows had MATLAB style indexing and not python. As a
  result the windows did not peak at the right place.
- Similarly, the window rounding was assymetric.
- Causes slight reduction in the reconstructed error in the toy data
  examples.
- `b` amplitudes now fall along PSD of the input data, fixing a
  lingering issue in the mrCOSTS fit.

Author: klapo

pydmd/costs.py

@@ -325,13 +325,11 @@ def calculate_lv_kern(window_length, corner_sharpness=None):
 
         lv_kern = (
             np.tanh(
-                corner_sharpness
-                * np.arange(1, window_length + 1)
-                / window_length
+                corner_sharpness * np.arange(0, window_length) / window_length
             )
             - np.tanh(
                 corner_sharpness
-                * (np.arange(1, window_length + 1) - window_length)
+                * (np.arange(0, window_length) - window_length - 1)
                 / window_length
             )
             - 1
@@ -354,7 +352,7 @@ def build_kern(window_length):
         """
         recon_filter_sd = window_length / 8
         recon_filter = np.exp(
-            -((np.arange(window_length) - (window_length + 1) / 2) ** 2)
+            -((np.arange(window_length) - (window_length - 1) / 2) ** 2)
             / recon_filter_sd**2
         )
         return recon_filter

pydmd/mrcosts.py

@@ -1041,14 +1041,8 @@ def global_scale_reconstruction(
 
             # Convolve each windowed reconstruction with a gaussian filter.
             # Std dev of gaussian filter
-            recon_filter_sd = mrd.window_length / 8
-            recon_filter = np.exp(
-                -(
-                    (np.arange(mrd.window_length) - (mrd.window_length + 1) / 2)
-                    ** 2
-                )
-                / recon_filter_sd**2
-            )
+            recon_filter = mrd.build_kern(mrd.window_length)
+
             omega_classes = omega_classes_list[n_mrd]
 
             if mrd.svd_rank < np.max(self._svd_rank_array):

tests/test_mrcosts.py

@@ -149,8 +149,8 @@ def rhs_UFD(t, y, eta, epsilon, tau):
 
 # Define the expected error in the reconstructions.
 expected_global_error = 0.053
-expected_lf_error = 0.12
-expected_hf_error = 0.19
+expected_lf_error = 0.10
+expected_hf_error = 0.17
 expected_transient_error = 0.3
 
 # Fit mrCOSTS for testing
@@ -360,7 +360,7 @@ def test_plot_local_time_series():
         _ = mrc.plot_local_time_series(0, 0, data=data.T)
 
 
-def tear_down():
+def test_tear_down():
     """Remove the files generated in `test_netcdf`"""
     file_list = glob.glob("*tests*.nc")
     for f in file_list: