Merge pull request PyDMD#532 from klapo/mrcosts-indexing-cleanup

Fixed bug in building convolution windows for mrCOSTS

Author: mtezzele

pydmd/costs.py

@@ -326,12 +326,11 @@ def calculate_lv_kern(window_length, corner_sharpness=None):
         lv_kern = (
             np.tanh(
                 corner_sharpness
-                * np.arange(1, window_length + 1)
-                / window_length
+                * np.linspace(0, 1, window_length, endpoint=False)
             )
             - 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 +353,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/meta.py

@@ -12,7 +12,9 @@
 
 __project__ = "PyDMD"
 __title__ = "pydmd"
-__author__ = "Nicola Demo, Marco Tezzele, Francesco Andreuzzi, Sara Ichinaga, Karl Lapo"
+__author__ = (
+    "Nicola Demo, Marco Tezzele, Francesco Andreuzzi, Sara Ichinaga, Karl Lapo"
+)
 __copyright__ = "Copyright 2017-2024, PyDMD contributors"
 __license__ = "MIT"
 __version__ = "1.0.0"

pydmd/mrcosts.py

@@ -560,7 +560,7 @@ def from_netcdf(self, file_list):
         self._n_data_vars = n_data_vars
         self._n_time_steps = n_time_steps
 
-    def to_netcdf(self, filename):
+    def to_netcdf(self, filename, filepath="."):
         """
         Save the mrCoSTS fit to file in netcdf format.
 
@@ -569,16 +569,22 @@ def to_netcdf(self, filename):
 
         :param filename: Common name shared by each file.
         :type filename: str
+        :param filepath: Path to save the results. Default is the current
+        directory.
+        :type filename: str
+
         """
         for c in self._costs_array:
+            fname = ".".join(
+                (
+                    filename,
+                    f"window={c.window_length:}",
+                    "nc",
+                )
+            )
+            fpath = os.path.join(filepath, fname)
             c.to_xarray().to_netcdf(
-                ".".join(
-                    (
-                        filename,
-                        f"window={c.window_length:}",
-                        "nc",
-                    )
-                ),
+                fpath,
                 engine="h5netcdf",
                 invalid_netcdf=True,
             )
@@ -1041,14 +1047,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
@@ -286,13 +286,14 @@ def test_omega_transforms():
         )
 
 
-def test_netcdf():
+def test_netcdf(tmp_path):
     """
     Test the round trip conversion of the mrCOSTS object to file in
     netcdf format and back to mrCOSTS.
     """
-    mrc.to_netcdf("tests")
-    file_list = glob.glob("*tests*.nc")
+    # Move the I/O tests to the temporary test directory.
+    mrc.to_netcdf("tests", filepath=tmp_path)
+    file_list = glob.glob(os.path.join(tmp_path, "*tests*.nc"))
     mrc_from_file = mrCOSTS()
     mrc_from_file.from_netcdf(file_list)
 
@@ -358,10 +359,3 @@ def test_plot_local_time_series():
 
     with raises(ValueError):
         _ = mrc.plot_local_time_series(0, 0, data=data.T)
-
-
-def tear_down():
-    """Remove the files generated in `test_netcdf`"""
-    file_list = glob.glob("*tests*.nc")
-    for f in file_list:
-        os.remove(f)