Fix the row order of H4 parameter files in MKZ_Param_Multi_Layer

CHANGELOG.md

@@ -2,10 +2,36 @@
 
 ## [Unreleased]
 
+- Changed
+  - `save_txt()` of `HH_Param_Multi_Layer` and `MKZ_Param_Multi_Layer` now
+    writes 6 significant digits by default (was 5), the same as the files that
+    `get_all_HH_x_params()` and `get_all_H4_x_params()` save
+  - `helper_mkz_model.deserialize_array_to_params()` now raises a `ValueError`
+    with a message when the array does not have 4 values, instead of a bare
+    `AssertionError` (which `python -O` skips)
+  - `MKZ_Param_Multi_Layer.serialize_to_2D_array()` now returns the layout of
+    the H4 parameter files that the class reads (rows gamma_ref, 0, s, beta),
+    instead of the rows gamma_ref, s, beta, Gmax. Code that reads rows of its
+    result by position needs to change. Like `save_txt()`, it warns when it
+    drops a Gmax other than 1.0 (#79)
 - Fixed
   - `Damping_Calibration.get_H4_x_param()` returning HH parameters
     (`HH_Param_Multi_Layer`) instead of H4 (MKZ) parameters
     (`MKZ_Param_Multi_Layer`) (#60)
+  - `MKZ_Param_Multi_Layer` reading H4 parameter files (and 2D arrays) in the
+    wrong row order, which made `s` 0 and the G/Gmax and damping curves flat.
+    Its `save_txt()` and `get_all_H4_x_params(save_txt=True)` now also write
+    the file format (gamma_ref, 0, s, beta). This format has no Gmax, so
+    `save_txt()` warns (once, naming the layers) when it drops a Gmax other
+    than 1.0. A file or 2D array that is not valid H4 parameters now raises a
+    `ValueError` that names the layer: its 2nd row ("b") is not 0, it has NaN
+    (for example, from a trailing delimiter), or it does not have 4 rows. H4
+    files written by 0.7.1 or earlier have the rows gamma_ref, s, beta, Gmax,
+    so they raise this error. To convert one, delete its 4th row and insert a
+    row of zeros as its 2nd row, e.g.,
+    `np.insert(np.genfromtxt(old_file, ndmin=2)[:3], 1, 0, axis=0)` (#79)
+  - `HH_Param_Multi_Layer` and `MKZ_Param_Multi_Layer` failing to load a
+    one-layer parameter file, including one written by `save_txt()` (#86)
 
 ## [0.7.1] - 2026-10-11
 

PySeismoSoil/class_curves.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import functools
 from pathlib import Path
 from typing import TYPE_CHECKING, Any, Literal
 
@@ -1186,9 +1187,12 @@ def get_all_H4_x_params(
             Number of CPU cores to use. If None, all cores are used. No effects
             if the parallelization options are set to ``False``.
         save_txt : bool, default=False
-            Whether to save the results as a "HH_x_STATION_NAME.txt" file.
+            Whether to save the results as a "H4_x_STATION_NAME.txt" file. The
+            file has one column per layer, with the rows gamma_ref, 0, s, beta
+            (the format that ``MKZ_Param_Multi_Layer`` reads). Gmax is not
+            saved.
         txt_filename : str | None, default=None
-            File name of the text file to save HH parameters. If the object is
+            File name of the text file to save H4 parameters. If the object is
             created via a "curve" text file, then ``txt_filename`` can be
             ``None`` and the output filename will be determined automatically.
         sep : str | None, default=None
@@ -1239,7 +1243,9 @@ def get_all_H4_x_params(
             save_txt=save_txt,
             txt_filename=txt_filename,
             sep=sep,
-            func_serialize=mkz.serialize_params_to_array,
+            func_serialize=functools.partial(
+                mkz.serialize_params_to_array, to_files=True
+            ),
         )
 
         return MKZ_Param_Multi_Layer(params)

PySeismoSoil/class_parameters.py

@@ -2,6 +2,7 @@
 
 import collections
 import json
+import warnings
 from collections.abc import Callable
 from typing import TYPE_CHECKING, Any, TypeAlias
 
@@ -506,23 +507,43 @@ def serialize_to_2D_array(self) -> np.ndarray:
         """
         Serialize the parameter data to a 2D numpy array.
 
+        The array has the layout of the parameter files that the subclass reads
+        (such as "HH_x_SITE_NAME.txt" or "H4_x_SITE_NAME.txt"; see the
+        docstring of the subclass), so ``save_txt()`` writes files that the
+        subclass can read back.
+
         Returns
         -------
         param_2D_array : np.ndarray
             A 2D numpy array whose columns are parameters of each layer.
         """
-        output = []
-        for param_single_layer in self.param_list:
-            param_array = param_single_layer.serialize()
-            output.append(param_array)
+        output = [self._serialize_layer(param) for param in self.param_list]
+        return np.array(output).T
+
+    @staticmethod
+    def _serialize_layer(param: Parameter) -> np.ndarray:
+        """
+        Serialize the parameters of one layer into one column of the 2D array.
+
+        Subclasses whose parameter files have a different layout than
+        ``Parameter.serialize()`` override this method.
 
-        param_2D_array = np.array(output).T
-        return param_2D_array  # noqa: RET504
+        Parameters
+        ----------
+        param : Parameter
+            The parameters of one layer.
+
+        Returns
+        -------
+        np.ndarray
+            The column of the 2D array for this layer.
+        """
+        return param.serialize()
 
     def save_txt(
             self,
             filename: str,
-            precision: str = '%.5g',
+            precision: str = hlp.PARAM_FILE_PRECISION,
             sep: str = '\t',
             **kw_to_savetxt: dict[Any, Any],
     ) -> None:
@@ -533,8 +554,10 @@ def save_txt(
         ----------
         filename : str
             File name (including path) of the output file.
-        precision : str, default='%.5g'
-            Precision of the numbers to be saved.
+        precision : str, default=hlp.PARAM_FILE_PRECISION
+            Precision of the numbers to be saved. The default (``'%.6g'``) is
+            the same as in the files that the curve-fitting methods (such as
+            ``get_all_HH_x_params()``) save.
         sep : str, default='\t'
             Delimiter identifier.
         **kw_to_savetxt : dict[Any, Any]
@@ -610,7 +633,8 @@ def __init__(
     ) -> None:
         if isinstance(filename_or_data, str):  # file name
             self._filename = filename_or_data
-            params = np.genfromtxt(filename_or_data, delimiter=sep)
+            # ``ndmin=2`` keeps a one-layer file (one column) 2D
+            params = np.genfromtxt(filename_or_data, delimiter=sep, ndmin=2)
             list_of_param_array = hlp.extract_from_param_format(params)
             list_of_param = [
                 hh.deserialize_array_to_params(_) for _ in list_of_param_array
@@ -658,20 +682,26 @@ class MKZ_Param_Multi_Layer(Param_Multi_Layer):
     Parameters
     ----------
     filename_or_data : MKZ_Param_Input
-        A file name of a validly formatted "parameter file", i.e., having the
+        A file name of a validly formatted "parameter file" (such as
+        "H4_G_SITE_NAME.txt" or "H4_x_SITE_NAME.txt"), i.e., having the
         following format:
             +----------------+-----------------+-----------------+-----+
             |  param_layer_1 |  param_layer_2  |  param_layer_3  | ... |
             +================+=================+=================+=====+
-            |      1.1       |      2.2        |      3.3        | ... |
+            |   gamma_ref    |   gamma_ref     |   gamma_ref     | ... |
             +----------------+-----------------+-----------------+-----+
-            |      1.2       |      2.3        |      3.4        | ... |
+            |       0        |       0         |       0         | ... |
             +----------------+-----------------+-----------------+-----+
-            |      ...       |      ...        |      ...        | ... |
+            |       s        |       s         |       s         | ... |
+            +----------------+-----------------+-----------------+-----+
+            |      beta      |      beta       |      beta       | ... |
             +----------------+-----------------+-----------------+-----+
 
         or a 2D numpy array containing the data of the format above, or a list
         containing MKZ parameter data (dictionaries or ``MKZ_Param`` objects).
+        The 2nd row ("b") is always 0, for historical reasons. This format has
+        no Gmax, so when reading a file or a 2D array, Gmax of each layer is
+        set to 1.0 (Gmax does not affect the G/Gmax or damping curves).
     sep : str, default='\t'
         Delimiter of the file to be imported. If ``filename_or_data`` is not a
         file name, ``sep`` has no effect.
@@ -687,6 +717,15 @@ class MKZ_Param_Multi_Layer(Param_Multi_Layer):
     ------
     TypeError
         When then type of ``filename_or_data`` is not valid
+    ValueError
+        When a column of the file or the 2D array is not valid H4 parameters:
+        it does not have 4 values, it has NaN (for example, from an empty cell
+        or a trailing delimiter in the file), or its 2nd value ("b") is not 0.
+        H4 files written by PySeismoSoil 0.7.1 or earlier (with
+        ``MKZ_Param_Multi_Layer.save_txt()`` or
+        ``get_all_H4_x_params(save_txt=True)``) have the rows gamma_ref, s,
+        beta, Gmax instead, so they raise this error. The error message says
+        how to convert them.
     """
 
     param_list: list[MKZ_Param]
@@ -700,24 +739,28 @@ def __init__(
     ) -> None:
         if isinstance(filename_or_data, str):  # file name
             self._filename = filename_or_data
-            params = np.genfromtxt(filename_or_data, delimiter=sep)
-            list_of_param_array = hlp.extract_from_param_format(params)
-            list_of_param = [
-                mkz.deserialize_array_to_params(_) for _ in list_of_param_array
-            ]
-        elif isinstance(filename_or_data, np.ndarray):
-            hlp.assert_2D_numpy_array(
-                filename_or_data, name='`filename_or_data`'
-            )
-            list_of_param_array = hlp.extract_from_param_format(
-                filename_or_data
-            )
-            list_of_param = [
-                mkz.deserialize_array_to_params(_) for _ in list_of_param_array
-            ]
-        elif isinstance(filename_or_data, list):
+            # ``ndmin=2`` keeps a one-layer file (one column) 2D
+            data = np.genfromtxt(filename_or_data, delimiter=sep, ndmin=2)
+        else:
+            data = filename_or_data
+
+        if isinstance(data, np.ndarray):
+            hlp.assert_2D_numpy_array(data, name='`filename_or_data`')
+            list_of_param = []
+            list_of_param_array = hlp.extract_from_param_format(data)
+            for j, param_array in enumerate(list_of_param_array):
+                try:
+                    param = mkz.deserialize_array_to_params(
+                        param_array, from_files=True
+                    )
+                except ValueError as e:
+                    msg = f'Layer #{j} (0-based): {e}'
+                    raise ValueError(msg) from e
+
+                list_of_param.append(param)
+        elif isinstance(data, list):
             self._filename = None
-            list_of_param = filename_or_data
+            list_of_param = data
         else:
             raise TypeError('Unrecognized type for ``filename_or_data``.')
 
@@ -727,3 +770,61 @@ def __init__(
             list_of_param,
             element_class=MKZ_Param,
         )
+
+    def serialize_to_2D_array(self) -> np.ndarray:
+        """
+        Serialize the parameter data to a 2D numpy array in the H4 file layout.
+
+        Each column is one layer, with the rows gamma_ref, 0, s, beta (see the
+        class docstring). This layout has no Gmax. If the Gmax of any layer is
+        not 1.0, one ``UserWarning`` names those layers (or says that it is
+        every layer).
+
+        Returns
+        -------
+        param_2D_array : np.ndarray
+            A 2D numpy array whose columns are parameters of each layer.
+        """
+        layers = [
+            j
+            for j, param in enumerate(self.param_list)
+            if param['Gmax'] != mkz.H4_FILE_GMAX
+        ]
+        if layers:
+            if len(layers) == len(self.param_list):
+                which_layers = 'Gmax is not 1.0 in every layer.'
+            else:
+                which_layers = (
+                    f'Layers whose Gmax is not 1.0 (0-based): {layers}.'
+                )
+
+            warnings.warn(
+                f'{mkz.H4_FILE_GMAX_WARNING} {which_layers}',
+                UserWarning,
+                stacklevel=hlp._find_stack_level(),  # noqa: SLF001
+            )
+
+        return super().serialize_to_2D_array()
+
+    @staticmethod
+    def _serialize_layer(param: Parameter) -> np.ndarray:
+        """
+        Serialize the parameters of one layer into the H4 file layout.
+
+        The layout is gamma_ref, 0, s, beta. Gmax is not included, so a file
+        written by ``save_txt()`` loads back with Gmax = 1.0 for each layer.
+        (``serialize_to_2D_array()`` has already warned about the layers whose
+        Gmax is not 1.0.)
+
+        Parameters
+        ----------
+        param : Parameter
+            The MKZ parameters of one layer.
+
+        Returns
+        -------
+        np.ndarray
+            The column of the 2D array for this layer.
+        """
+        param_without_Gmax = {**param.data, 'Gmax': mkz.H4_FILE_GMAX}
+        return mkz.serialize_params_to_array(param_without_Gmax, to_files=True)

PySeismoSoil/helper_generic.py

@@ -2,8 +2,11 @@
 
 from __future__ import annotations
 
+import inspect
+import os
 import platform
 from datetime import datetime
+from pathlib import Path
 from typing import TYPE_CHECKING, Any, Literal
 
 import matplotlib.pyplot as plt
@@ -20,6 +23,10 @@
 NDIM_2D_ARRAY = 2
 NUM_COLUMNS_TWO_COL_DATA = 2
 
+# Precision of the numbers in the parameter files (such as "HH_x_SITE.txt" or
+# "H4_x_SITE.txt") that PySeismoSoil writes
+PARAM_FILE_PRECISION = '%.6g'
+
 # A Vs profile has either 2 columns (thickness, Vs) or 5 columns (thickness,
 # Vs, damping, density, material number)
 NUM_COLUMNS_THICKNESS_AND_VS = 2
@@ -36,6 +43,34 @@
 CHECK_STATUS_NEGATIVE = -3
 
 
+def _find_stack_level() -> int:
+    """
+    Find the ``stacklevel`` for ``warnings.warn()`` that points at user code.
+
+    A warning issued deep inside PySeismoSoil (for example, from
+    ``serialize_to_2D_array()`` when it is called by ``save_txt()``) would
+    otherwise point at a line of this library. This function counts the frames
+    inside the PySeismoSoil package, starting from the function that calls it,
+    so that the warning points at the first line outside of it.
+
+    Returns
+    -------
+    int
+        The ``stacklevel`` to pass to ``warnings.warn()``.
+    """
+    package_dir = f'{Path(__file__).parent}{os.sep}'
+    frame = inspect.currentframe()
+    frame = frame.f_back if frame is not None else None  # the caller
+    level = 1
+    while frame is not None and frame.f_code.co_filename.startswith(
+        package_dir
+    ):
+        frame = frame.f_back
+        level += 1
+
+    return level
+
+
 def detect_OS() -> str:
     """
     Check which operating system is currently running.