Fix the row order of H4 parameter files in MKZ_Param_Multi_Layer

[jsh9]

Fixes #79. Fixes #86. One of the four P0 fixes from #114.

Problem

H4 parameter files (H4_G_*.txt, H4_x_*.txt) store each layer as gamma_ref, b (always 0), s, beta. MATLAB SeismoSoil writes this layout, and its H4 kernel reads it (see "Where the layout comes from" below).

• Reading: MKZ_Param_Multi_Layer() read files and 2D arrays without from_files=True, so every layer got s = 0, beta = the file's s, and Gmax = the file's beta. With s = 0, G/Gmax is a constant (0.5296 for every strain in layer 1 of H4_G_IWTH04.txt), and damping is about 0.
• Writing: save_txt() and get_all_H4_x_params(save_txt=True) wrote gamma_ref, s, beta, Gmax instead of the file layout.

Reading and writing were wrong in the same way, so a file written by the library loaded back unchanged. That is why the existing round-trip tests passed.

Separately, HH_Param_Multi_Layer() and MKZ_Param_Multi_Layer() could not load a one-layer parameter file (#86): np.genfromtxt() returns a 1D array for a file with one column.

Changes

• PySeismoSoil/class_parameters.py:
  • MKZ_Param_Multi_Layer.__init__() reads files and 2D arrays with from_files=True. The class docstring said a 2D array has the same format as a file, but its table did not name the rows. It now does. The file branch loads the file into a 2D array and falls through to the 2D-array branch.
  • A file or 2D array that is not valid H4 parameters raises a ValueError whose message starts with the first bad layer (Layer #j (0-based): ...): its 2nd row ("b") is not 0, it has NaN (for example, from an empty cell or a trailing delimiter), or it does not have 4 rows (for example, an HH file). MATLAB's H4 fitters always write 0 in the 2nd row (see "Where the layout comes from"), so their files load. H4 files written by PySeismoSoil 0.7.1 or earlier have s there, so they raise this error instead of loading wrongly, and the message says how to convert them. The checks themselves are in mkz.deserialize_array_to_params() (below), so they also apply to direct callers.
  • HH_Param_Multi_Layer and MKZ_Param_Multi_Layer read files with np.genfromtxt(..., ndmin=2), so one-layer files load (HH_Param_Multi_Layer() and MKZ_Param_Multi_Layer() cannot load a one-layer parameter file, including one written by save_txt() #86).
  • Param_Multi_Layer.serialize_to_2D_array() calls a per-layer hook, _serialize_layer() (default: Parameter.serialize()). MKZ_Param_Multi_Layer overrides the hook to write the file layout (to_files=True), so save_txt() does too.
  • The file layout has no Gmax: reading sets Gmax = 1.0 and writing drops it. This is in the docstrings. Gmax does not affect the G/Gmax or damping curves. When the dropped Gmax is not 1.0, MKZ_Param_Multi_Layer.serialize_to_2D_array() (and so save_txt()) issues one UserWarning instead of one warning per layer. It ends with "Layers whose Gmax is not 1.0 (0-based): [0, 2]." or, when every layer is affected (for example, Gmax from a Vs profile), "Gmax is not 1.0 in every layer." It then calls the base method, which serializes each layer with Gmax = 1.0, so the per-layer helper does not warn again.
  • save_txt() now defaults to hlp.PARAM_FILE_PRECISION ('%.6g', was '%.5g'), the same as the files that the curve-fitting methods write.
• PySeismoSoil/helper_mkz_model.py:
  • deserialize_array_to_params() raises a ValueError with a message when the array does not have 4 values (it was a bare assert, which python -O skips). With from_files=True, it also raises a ValueError for NaN or infinite values, and for a "b" that is not within H4_FILE_B_TOLERANCE (1e-8) of 0.
  • serialize_params_to_array(to_files=True) issues a UserWarning when Gmax is not 1.0, because the file layout cannot store it. The warning text is shared with MKZ_Param_Multi_Layer as H4_FILE_GMAX_WARNING.
• PySeismoSoil/helper_generic.py:
  • New constant PARAM_FILE_PRECISION = '%.6g'.
  • New private helper _find_stack_level(). It counts the frames inside the package, so that a warning points at the first line outside it: the user's save_txt() or serialize_to_2D_array() call, not a line in class_parameters.py. Both Gmax warnings use it, with # noqa: SLF001 like the other cross-module private helpers (e.g., _process_fig_ax_objects()). It is private so that the API docs do not list it.
• PySeismoSoil/class_curves.py:
  • get_all_H4_x_params() serializes with functools.partial(mkz.serialize_params_to_array, to_files=True).
  • Its docstring said it saves an "HH_x" file of "HH parameters"; it now says H4 and describes the layout.
• PySeismoSoil/helper_site_response.py: fit_all_damping_curves() writes with hlp.PARAM_FILE_PRECISION, and its func_serialize docstring mentions to_files=True for H4 files.
• CHANGELOG.md: "Fixed" entries for MKZ_Param_Multi_Layer() reads H4 parameter files in the wrong row order, so s is 0 and the G/Gmax curves are flat #79 (including that the Gmax warning is issued once and names the layers, and a one-line recipe to convert old H4 files) and HH_Param_Multi_Layer() and MKZ_Param_Multi_Layer() cannot load a one-layer parameter file, including one written by save_txt() #86, and "Changed" entries for the save_txt() precision, the ValueError in deserialize_array_to_params(), and the new row layout of MKZ_Param_Multi_Layer.serialize_to_2D_array() (code that reads its rows by position needs to change).
• examples/Demo_05_Parameters.ipynb: it said MKZ_Param_Multi_Layer behaves identically to HH_Param_Multi_Layer. It now says that its files, serialize_to_2D_array(), and save_txt() use the H4 file layout, which has no Gmax.
• The example notebooks were re-run (ae4e873) with Python 3.12, the version recorded on main, as CONTRIBUTING.md requires. No notebook loads H4 files, builds MKZ_Param_Multi_Layer from a file or array, or calls save_txt(), so no printed results changed. Apart from the Demo_05 note, only timings, folder names, object addresses, progress bars, and image bytes changed. The later commits (6bf6beb, 4dfa939, 3da325f) change no notebook output: no notebook triggers the Gmax warning.

Where the layout comes from

PySeismoSoil. The docstrings of serialize_params_to_array() and deserialize_array_to_params() in helper_mkz_model.py describe both orders. The first says "the text files recognizable by MATLAB and Fortran functions have the convention of "gamma_ref, 0.0, s, beta"". The library's own MKZ fit already uses this layout. fit_MKZ(), a port of MATLAB's fitMKZ.m with the same bounds, returns one row per layer:

    param = np.column_stack((ref_strain, np.zeros(n_ma), s_value, beta))

The HH calibration (helper_hh_calibration.py) reads it that way:

        param_MKZ, _ = mkz.fit_MKZ(curves, show_fig=show_fig)
        gamma_ref = param_MKZ[:, 0]
        s = param_MKZ[:, 2]
        beta = param_MKZ[:, 3]

The git history agrees. 4b80e2a (2019-02-15) added both orders: gamma_ref, s, beta, Gmax for in-memory serialization (kept by MKZ_Param.serialize(), which this PR does not change) and gamma_ref, 0, s, beta for files. The first get_all_H4_x_params() wrote its file with to_files=True; the refactor in 862f996 (2019-03-02) dropped that flag. 2838e75 (2019-04-11) added H4_G_IWTH04.txt (2nd row all 0) as the test input of MKZ_Param_Multi_Layer, and 72ebd9e (same day) added a test that serialize_to_2D_array() equals that file. That test passed only because reading and writing used the same wrong order.

MATLAB SeismoSoil (Caltech-geoquake/SeismoSoil). tests/files/H4_G_IWTH04.txt and H4_x_IWTH04.txt are byte-for-byte copies of its Sample_input_files/. The GUI writes the fitted matrices to the H4 files as they are (SeismoSoil_Preparation_H2_H4_Curve_Fitting.m):

    H4G = fitMKZ(curve,handles.metricdata.showfig);
    ...
    H4x = gaH4x(curve,handles.metricdata.showfig);
    ...
    dlmwrite(fullfile(curve_dir_name,H4G_filename),H4G,'delimiter','\t');
    dlmwrite(fullfile(curve_dir_name,H4x_filename),H4x,'delimiter','\t');

runNLH4FromGUI.m writes the imported matrices for the kernel as they are, then runs NLH4.exe or nlH4.m:

dlmwrite('H4_G.dat', para_G);
dlmwrite('H4_x.dat', para_x);

MATLAB source of the functions that SeismoSoil ships as compiled .p files (not public). fitMKZ.m makes the H4_G parameters:

% Note: para = [gamma_ref,  0, s, beta]
%                ^ unit: 1  ^ This is always zero...
%                            (Just for compatibility with legacy codes)
...
para=[ref_strain; zeros(1,n_ma); s_value; beta];
...
opts.Lower = [0.2, 0, 0.6];  % define lower bounds for: [beta, gamma_ref, s]
opts.Upper = [1.8, 0.1, 0.999];  % upper bounds for: [beta, gamma_ref, s]

gaH4x.m makes the H4_x parameters. Its model has a strain-dependent reference strain, x(1)+gamma*x(2), but it uses only the first column of LB/UB, which fixes b (row 2) at 0. gaH4G.m has the same bounds and the same output line.

LB=[0.0005, 0.000001, 0.000001, 0.000001; ...
    0,      0.005,     0.005,      0.005; ...
    0.6, 0.7, 0.7, 0.7; ...
    0.2 0.5 0.5 0.5 ];
UB=[0.005, 0.01, 0.01, 0.01; ...
    0,     1.5,    1.5,  1.5;  ...
    0.999, 1.3, 1.3, 1.3; ...
    1.8 5 5 5];
...
    [ref_strain(k),const_b(k),s_value(k),beta(k)] = ...
        get_parameter_H2_x(strn_G_m(:,k), G_vec_m(:,k), strn_x_m(:,k), ...
                           x_vec_m(:,k), LB(:,1), UB(:,1), nr_iter, verbosity_level);
...
para=[ref_strain; const_b; s_value; beta];
...
function [ratio,damping]=my_H2(gamma,x)
strain_ref=x(1)+gamma*x(2);
ratio=1/(1+x(4)*(gamma/strain_ref)^x(3));

nlH4.m, the MATLAB version of the H4 kernel, reads rows 1, 3 and 4 as gamma_ref, s and beta, and does the same for H4_x.dat. It reads row 2 into b_val, but no calculation uses it. Gmax comes from the Vs profile, not from the file.

    G(i) = r(i) * v(i) * v(i);
...
para_H2 = importdata('H4_G.dat');
...
    a_val(k) = para_H2(1,ma(k));
    b_val(k) = para_H2(2,ma(k));
    s_val(k) = para_H2(3,ma(k));
    beta(k) = para_H2(4,ma(k));
...
        Gmax(1,lay) = G(i);
...
    gamma_ref = a_val_sub(G_x,mat);
    fbb = Gmax(G_x,mat)*xx/(1.0+beta_sub(G_x,mat)*(abs(xx/gamma_ref))...
        ^s_val_sub(G_x,mat));

The sample files fit these bounds. Row 3 is 0.886–0.889 (H4_G) and 0.610–0.686 (H4_x), inside the bounds for s (0.6–0.999). Row 4 is 0.47–1.79 and 0.48–1.41, inside the bounds for beta (0.2–1.8). main reads these files with s = 0, below the lower bound.

Tests

These fail on the old library code (main) and pass with the fix:

• test_mkz_param_multi_layer__test_contents_of_list_elements: the old benchmark [0.000856, 0, 0.88832, 1.7492] encoded s = 0. The test now expects layer 1 of H4_G_IWTH04.txt read in its file layout: gamma_ref=0.00085636, s=0.88832, beta=1.7492, Gmax=1.0.
• test_mkz_param_multi_layer__save_txt_and_load_back (new): serialize_to_2D_array() and the save_txt() file have the rows gamma_ref, 0, s, beta, both warn that a Gmax other than 1.0 is not saved, and both load back with the same parameters (Gmax = 1.0).
• test_mkz_param_multi_layer__one_Gmax_warning_at_the_caller[serialize_to_2D_array|save_txt] (new): with Gmax other than 1.0 in layers 0 and 2 of 3, each method issues exactly one warning, which ends with "Layers whose Gmax is not 1.0 (0-based): [0, 2]." and points at the test file (the caller).
• test_mkz_param_multi_layer__Gmax_warning_if_every_layer_is_affected[1|3] (new): when every layer (1 of 1, or 3 of 3) has a Gmax other than 1.0, there is one warning, which ends with "Gmax is not 1.0 in every layer." instead of listing the layers.
• test_construct_curves[from_H4_G_parameters] / [from_H4_x_parameters]: these checked only shapes. They now also compare the values of one layer with the curve from that layer's column read with from_files=True. The HH cases get the same check.
• test_mkz_param_multi_layer__invalid_H4_parameters_raise_error[file|2D_array × old_format|one_layer_with_nonzero_b|NaN|too_few_rows|HH_parameters] (new): each raises the ValueError, and the message names the bad layer.
• test_mkz_param_multi_layer__trailing_delimiter_raises_error (new): a file with a trailing tab on each line raises the NaN error for the extra column, not the "b" error.
• test_HH_MKZ_param_multi_layer__one_layer_file_loads_back[HH|MKZ] (new, HH_Param_Multi_Layer() and MKZ_Param_Multi_Layer() cannot load a one-layer parameter file, including one written by save_txt() #86): a one-layer object saved with save_txt() loads back from the file.
• tests/test_class_curves.py::test_get_all_H4_x_params__saved_file_loads_back (new): the saved file has the file layout, loads back via MKZ_Param_Multi_Layer(), and is byte-for-byte the same as the file save_txt() writes (same layout and precision).
• tests/test_helper_mkz_model.py:
  • test_deserialize_array_to_params__failure[incorrect_number_of_parameters] now expects a ValueError with a message instead of AssertionError.
  • test_serialize_params_to_array__to_files_warns_if_Gmax_is_not_1 (new): it also checks that the warning points at the caller.
  • test_deserialize_array_to_params__from_files_failure[old_layout|NaN|incorrect_number_of_parameters] (new).

These check that nothing else changed, and pass on both: test_mkz_param_multi_layer__save_txt_without_warning_if_Gmax_is_1, test_mkz_param_multi_layer__b_close_to_0_is_accepted, test_serialize_params_to_array__to_files, and test_deserialize_array_to_params__from_files[b_is_0|b_is_close_to_0].

The single-layer H4 and HH tests now build their parameters with the _H4_param_from_column() / _HH_param_from_column() helpers.

Testing

Run locally on 4dfa939, the last commit that changes code (3da325f changes only CHANGELOG.md):

• Python 3.13, python -m pytest tests: 269 passed. Against the library code of 6bf6beb (the head before 4dfa939), the 2 cases of test_mkz_param_multi_layer__one_Gmax_warning_at_the_caller and the 2 cases of test_mkz_param_multi_layer__Gmax_warning_if_every_layer_is_affected fail.
• Python 3.12: tests/test_class_parameters.py, tests/test_helper_mkz_model.py, and tests/test_class_curves.py pass (75 passed).
• pydoclint PySeismoSoil: pass. pre-commit run --all-files: pass, also on 3da325f.
• The CHANGELOG recipe converts an H4 file written by main's save_txt() into the new layout, and MKZ_Param_Multi_Layer() then loads the original parameters. Without the conversion, the file raises the ValueError, whose message explains the conversion.

Notes for review

• Behavior change for files written before this fix: H4 files written by get_all_H4_x_params(save_txt=True) or MKZ_Param_Multi_Layer.save_txt() in 0.7.1 or earlier are in the old layout (gamma_ref, s, beta, Gmax). They now raise a ValueError instead of loading wrongly. The message says how to convert them, and the CHANGELOG has a one-line recipe. The Gmax they stored is 1.0 for get_all_H4_x_params() (the H4_x fit always sets it), so converting them loses nothing. Files in the H4 layout, such as the ones in tests/files/ and those made by MATLAB SeismoSoil, now load correctly.
• save_txt() precision: the default is now 6 significant digits for HH_Param_Multi_Layer too.
• serialize_params_to_array() and deserialize_array_to_params() check their input length with assert, so the check disappears under python -O #89 is only partly covered: mkz.deserialize_array_to_params() no longer uses assert for its length check. mkz.serialize_params_to_array() and both HH functions still do.
• Not changed here:
  • Nonlinear_Simulation accepts MKZ_Param_Multi_Layer for G_param/xi_param according to its docstring. But it writes serialize_to_2D_array() into HH_G.dat/HH_x.dat and runs the HH kernel (PySeismoSoil ships only NLHH), which reads 9 HH parameters per layer. So MKZ parameters do not work there with either layout: the kernel fails to read HH_G.dat, and run() then fails with an unrelated ValueError (Nonlinear_Simulation accepts MKZ_Param_Multi_Layer, as documented, but run() crashes because the NLHH kernel reads only HH parameters #124).
  • HH_Param_Multi_Layer() does not reject NaN parameters (from a trailing delimiter or an empty cell), unlike MKZ_Param_Multi_Layer() here (HH_Param_Multi_Layer() accepts NaN parameters, so a trailing delimiter adds a layer of NaN and an empty cell gives NaN curves, without an error #125).
  • The damping fitters still change the curve data in place (fit_HH_x_single_layer() and fit_H4_x_single_layer() subtract the initial damping in place, changing the Damping_Curve's data and the caller's array #84).
• Fix Damping_Calibration.get_H4_x_param() returning HH parameters #115 (Damping_Calibration.get_H4_x_param() returns HH parameters instead of H4 (MKZ) parameters #60) is merged, and main is merged into this branch. Make the filters zero-phase again and use padlen #117 and Apply the transfer function at the right FFT bins in amplify_motion() #118 change other library files. All three re-run the notebooks and add an entry to the "Fixed" group in CHANGELOG.md (this PR also adds a "Changed" group), so after one of them merges, the others need main merged in again (keep both entries, then re-run the notebooks).

🤖 Generated with Claude Code

https://claude.ai/code/session_01W7fK8T2pPLf3tPCnMCKL1T
https://claude.ai/code/session_01FSqLqxJej1HeeayXnFVrvV
https://claude.ai/code/session_018XWkqdYHdQvARgx5RDCmcC