Refactor tabular interpolation

.github/workflows/ci.yml

@@ -13,7 +13,7 @@ jobs:
       - uses: actions/checkout@v3
       - uses: actions/setup-python@v4
         with:
-          python-version: "3.10"
+          python-version: "3.11"
           cache: pip
       - name: Run pre-commit
         uses: pre-commit/action@v3.0.0
@@ -26,7 +26,7 @@ jobs:
           - ubuntu-latest
           # - windows-latest Reason: get it running on ubuntu-latest first
           # - macOS-latest   Reason: get it running on ubuntu-latest first
-        python-version: ["3.10",]
+        python-version: ["3.11",]
     name: Run tests
     runs-on: ${{ matrix.os }}
     steps:

.gitignore

@@ -11,6 +11,7 @@
 __pycache__/
 *.py[cod]
 *$py.class
+/bin
 
 # C extensions
 *.so

README.md

@@ -53,7 +53,7 @@ The corresponding implementation is available in our [AdaptiveThermoMechROM](htt
 
 ### Online phase: Usage of the thermo-mechanical NTFA in simulations on the macroscale
 
-1. Load the tabular data for the NTFA matrices $`\underline{\underline{A}}(\theta_j)`$, $`\underline{\underline{D}}(\theta_j)`$, $`\bar{\underline{\underline{C}}}(\theta_j)`$, $`\underline{\tau}\_{\mathrm{\theta}}(\theta_j)`$, and $`\underline{\tau}_{\mathsf{p}}(\theta_j)`$ that are generated in the offline phase based on direct numerical simulations on the microscale.
+1. Load the tabular data for the NTFA matrices $`\underline{\underline{A}}(\theta_j)`$, $`\underline{\underline{D}}(\theta_j)`$, $`\bar{\underline{\underline{C}}}(\theta_j)`$, $`\underline{\tau}_{\mathrm{\theta}}(\theta_j)`$, and $`\underline{\tau}_{\mathsf{p}}(\theta_j)`$ that are generated in the offline phase based on direct numerical simulations on the microscale.
 Optionally truncate the NTFA modes $`N_{\mathrm{modes}}`$ to be used.
 
 2. Perform a linear interpolation to determine the NTFA matrices at the current model temperature based on the tabular data.
@@ -86,7 +86,11 @@ pip install -e .
 
 ### Requirements
 
-TODO: update
+TODO: update dependencies
+
+TODO: upload dataset to Darus
+
+TODO: provide functionality for download from Darus
 
 - Python 3.9 or later
 - `numpy` and `h5py` (installed as part of the `thermontfa` PIP package)

docs/README.md

@@ -0,0 +1,19 @@
+# Building the ThermoNTFA Python documentation
+
+To build the documentation:
+
+1. Install ThermoNTFA (`thermontfa` pip package). It must be possible to import
+   the module ``thermontfa``.
+2. Run in this directory:
+
+       make html
+
+## Processing demo/example programs
+
+Python demo programs are written in Python, with MyST-flavoured
+Markdown syntax for comments.
+
+1. `jupytext` reads the Python demo code, then converts to and writes a
+   Markdown file.
+2. `myst_parser` allows Sphinx to process the Markdown file.
+

test/test_interpolation.py

@@ -3,36 +3,35 @@
 Test tabular interpolation
 """
 
-# %%
 import h5py
 import numpy as np
-
-# %%
 from thermontfa import TabularInterpolation
 
 
-# %%
 def test_tab_inter_init():
+    theta = np.linspace(0, 1, 10)
+    val = np.linspace(2, 5, 10).reshape((-1, 1))
+
+    tab_inter = TabularInterpolation(temps=theta, data=val)
+    assert tab_inter.temp_min == theta.min()
+    assert tab_inter.temp_max == theta.max()
+    assert np.allclose(tab_inter.temps.ravel(), theta.ravel())
+    assert np.allclose(tab_inter.data.ravel(), val.ravel())
+
+
+def test_tab_inter_from_h5():
     theta = np.arange(5)
     val = theta**2
     idx = np.random.permutation(np.arange(5))
-    theta = theta[idx]
-    val = val[idx]
-    val.reshape((5, 1))
-    F = h5py.File("test.h5", "w")
-    F.create_dataset("/theta", data=theta)
-    F.create_dataset("/data", data=val)
-    F.close()
-
-    tab_inter = TabularInterpolation()
-    tab_inter.init_h5("test.h5", "/theta", "/data")
-    print(tab_inter.t_min)
-    print(tab_inter.t_max)
-    print(tab_inter.t)
-    print(tab_inter.data)
+    theta_perm = theta[idx]
+    val_perm = val[idx]
+    val_perm = val_perm.reshape((-1, 1))
+    with h5py.File("test.h5", "w") as file:
+        file.create_dataset("/theta", data=theta_perm)
+        file.create_dataset("/data", data=val_perm)
 
-
-# %%
-a = np.linspace(2, 5, 10).reshape((-1, 1))
-t = np.linspace(0, 1, 10)
-inter = TabularInterpolation(a, t)
+    tab_inter = TabularInterpolation.from_h5(file_name="test.h5", dset_temps="/theta", dset_data="/data")
+    assert tab_inter.temp_min == theta.min()
+    assert tab_inter.temp_max == theta.max()
+    assert np.allclose(tab_inter.temps.ravel(), theta.ravel())
+    assert np.allclose(tab_inter.data.ravel(), val.ravel())

thermontfa/tabular_interpolation.py

@@ -19,125 +19,138 @@
 390740016 (EXC-2075); 406068690 (FR2702/8-1); 517847245 (FR2702/10-1).
 """
 
-from typing import Optional, Tuple
+from typing import Optional, Tuple, Self
 
 import h5py
 import numpy as np
 
 
 class TabularInterpolation:
     """
-    Tabular interpolation for the NTFA
+    Tabular interpolation for the thermo-mechanical NTFA
+
+    Performs a linear interpolation of the NTFA system matrices for a given temperature
+    given tabular data for sufficiently many temperature points.
+    It can be initialized for given data or based on a HDF5 file.
     """
 
-    t_min: float = 0.0
-    t_max: float = 0.0
+    temp_min: float = 0.0
+    temp_max: float = 0.0
     dim: Tuple[int, ...] = ()
-    t: np.ndarray = np.array([])
+    temps: np.ndarray = np.array([])
     data: np.ndarray = np.array([])
     const_extrapolate: bool = False
 
     def __init__(
-        self, data: Optional[np.ndarray] = None, param: Optional[np.ndarray] = None
-    ):
+        self, temps: np.ndarray = None, data: np.ndarray = None, const_extrapolate: bool = False
+    ) -> None:
         """
-        Initialize the tabular interpolator
-
-        :param data:
-        :param param:
+        Initialize the tabular interpolator for given `data` at prescribed temperatures `temps`.
+
+        :param temps: temperature points on which tabular data is available.
+            The shape of the `numpy` array is expected to be `(N_t,)`.
+        :type temps: np.ndarray
+        :param data: tabular data, e.g., a batch of NTFA system matrices with shape `(N_t, ...)`.
+        :type data: np.ndarray
+        :param const_extrapolate: If true, a constant extrapolation instead of a linear extrapolation is performed.
+            The default value is false.
+        :type const: bool
         """
-        if data is not None and param is not None:
-            idx = np.argsort(param)
+        if temps is None or data is None:
+            raise ValueError("The arguments `temps` and `data` are required.")
+
+        idx = np.argsort(temps)
+
+        self.data = data[idx, ...]
+        self.temps = temps[idx]
+        self.temp_min = self.temps[0]
+        self.temp_max = self.temps[-1]
+        self.dim = data[0].shape
+        self.const_extrapolate = const_extrapolate
 
-            self.data = data[idx, :]
-            self.t = param[idx]
-            self.t_min = self.t[0]
-            self.t_max = self.t[-1]
-            self.dim = data[0].shape
+        n = self.temps.shape[0]
 
-    def init_h5(
-        self,
+        if data.ndim == 1:
+            self.data = self.data.reshape((-1, 1))
+        assert (
+            self.temps.size == n
+        ), "Parameter must be an array of scalars (i.e. shape=[n] or [n, 1] or [n, 1, 1] or ...)"
+        assert (
+            self.data.shape[0] == n
+        ), f"Number of scalar parameters not matching dimension of available data ({n} vs. {self.data.shape[0]}."
+
+    @classmethod
+    def from_h5(
+        cls,
         file_name: str,
-        dset_param: str,
+        dset_temps: str,
         dset_data: str,
-        extrapolate: str = "linear",
-        transpose: Optional[Tuple[int, ...]] = None,
-    ):
+        transpose_dims: Optional[Tuple[int, ...]] = None,
+        const_extrapolate: bool = False,
+    ) -> Self:
         """
-        Initialize the tabular interpolator using data in a H5 file
-
-        :param file_name: path of the H5 file
+        Initialize the tabular interpolator based on tabular data stored in a HDF5 file (*.h5).
+
+        This is a factory method and returns a new instance of the `TabularInterpolation` class.
+        It is expected that the HDF5 file contains a data set with path `dset_temps` that contains
+        a list of the temperature points on which tabular data is available. The shape of this
+        dataset is expected to be `(N_t,)`.
+        Additionaly, the HDF5 file must contain a data set with path `dset_data` that contains the
+        tabular data, e.g., a batch of NTFA system matrices with shape `(N_t, ...)`.
+        The order of axes/dimensions of the data set with path `dset_data` can be changed by transposing
+        to the axis order given in `transpose_dims`.
+
+        :param file_name: path of the HDF5 file
         :type file_name: str
-        :param dset_param: path to the desired dataset in the H5 file
+        :param dset_param: path to the desired dataset in the HDF5 file
         :type dset_param: str
         :param dset_data:
         :type dset_param: str
-        :param extrapolate: "linear" or "constant"
-        :type extrapolate: str
-        :param transpose: axis order for transposition
-        :type transpose: Tuple[int, ...], optional
+        :param const_extrapolate: "linear" or "constant"
+        :type extrapolate: bool
+        :param transpose_dims: axis order for transposition
+        :type transpose_dims: Tuple[int, ...], optional
+        :return: new instance  of the `TabularInterpolation` class
+        :rtype: TabularInterpolation
         """
-        F = h5py.File(file_name, "r")
-        self.t = np.array(F[dset_param])
-        if transpose is None:
-            self.data = np.array(F[dset_data])
-        else:
-            self.data = np.array(F[dset_data]).transpose(transpose)
-        F.close()
-        self.dim = self.data.shape[1:]
-
-        if extrapolate == "linear":
-            self.const_extrapolate = False
-        else:
-            if extrapolate == "constant":
-                self.const_extrapolate = True
+        with h5py.File(file_name, "r") as file:
+            temps = np.array(file[dset_temps])
+
+            if transpose_dims is None:
+                data = np.array(file[dset_data])
             else:
-                raise ValueError(
-                    f'extrapolate can only take values "linear" or "constant" (received: "{extrapolate}"'
-                )
+                data = np.array(file[dset_data]).transpose(transpose_dims)
+
+        return cls(temps=temps, data=data, const_extrapolate=const_extrapolate)
 
-        if self.data.ndim == 1:
-            self.data = self.data.reshape((-1, 1))
-        n = self.t.shape[0]
-        assert (
-            self.t.size == n
-        ), "ERROR: parameter must be an array of scalars (i.e. shape=[n] or [n, 1] or [n, 1, 1] or ...)"
-        assert (
-            self.data.shape[0] == n
-        ), f"ERROR: number of scalar parameters not matching dimension of available data ({n} vs. {self.data.shape[0]}."
-        idx = np.argsort(self.t)
-        self.t = self.t[idx]
-        self.data = self.data[idx, :]
-        self.t_min = self.t[0]
-        self.t_max = self.t[-1]
-
-    def interpolate(self, t: float) -> np.ndarray:
+    def interpolate(self, temp: float) -> np.ndarray:
         """
-        Perform a linear interpolation at a given temperature t
+        Perform a linear interpolation based on the available tabular data at a given temperature `temp`
 
-        :param t: temperature point for interpolation
-        :type t: float
+        :param temp: temperature point for interpolation
+        :type temp: float
         :return: interpolated quantity
         :rtype: np.ndarray
         """
-        if t < self.t_min:
+        if temp < self.temp_min:
             if self.const_extrapolate:
-                return self.data[0, :]
+                return self.data[0, ...]
             else:
                 # linear extrapolation
-                alpha = (self.t_min - t) / (self.t[1] - self.t[0])
-                return self.data[0, :] - alpha * (self.data[1, :] - self.data[0, :])
+                alpha = (self.temp_min - temp) / (self.temp[1] - self.temp[0])
+                return self.data[0, ...] - alpha * (self.data[1, ...] - self.data[0, ...])
 
-        if t >= self.t_max:
+        if temp >= self.temp_max:
             if self.const_extrapolate:
-                return self.data[-1, :]
+                return self.data[-1, ...]
             else:
                 # linear extrapolation
-                alpha = (t - self.t_max) / (self.t[-1] - self.t[-2])
-                return self.data[-1, :] + alpha * (self.data[-1, :] - self.data[-2, :])
+                alpha = (temp - self.temp_max) / (self.temp[-1] - self.temp[-2])
+                return self.data[-1, ...] + alpha * (self.data[-1, ...] - self.data[-2, ...])
 
-        idx = np.searchsorted(self.t > t, 1) - 1
+        idx = np.searchsorted(self.temp > temp, 1) - 1
         t1 = self.t[idx]
         t2 = self.t[idx + 1]
-        alpha = (t - t1) / (t2 - t1)
-        return self.data[idx, :] + alpha * (self.data[idx + 1, :] - self.data[idx, :])
+        alpha = (temp - t1) / (t2 - t1)
+        interp_data = self.data[idx, ...] + alpha * (self.data[idx + 1, ...] - self.data[idx, ...])
+        return interp_data