Goodness-of-Fit Scoring, QOL, and Functional Updates (#22)

* Fix multiple curve plotting

Remove mandatory creation of new fig and ax in Multiple_Curves.plot() to allow for superclass Curve's plot() function to handle plotting according to documentation

* Revert "Fix multiple curve plotting"

This reverts commit f5e4873.

* Fix multiple curve plotting

Remove mandatory creation of new fig and ax in Multiple_Curves.plot() to make functionality consistent with documentation

* Modify rho calculation for edge surface layer case

Prevent unphysical rho calculations due to low overburden pressure for a thin surface layer

* Update filter kernel to use SOS representation

As per scipy.signal.butter, "If the transfer function form [b, a] is requested, numerical problems can occur since the conversion between roots and the polynomial coefficients is a numerically sensitive operation"

* Fix nonlinear simulation dummy curve curve issue

Even though the dummy curves are not simulation significant, since the ggmax and damping parameters are calibrated separately, allows for the proper curves to be created by their corresponding calibrated parameters to avoid triggering a ValueError related to unrealistic damping curves created from the parameters calibrated for the ggmax curves

* Batch nonlinear simulation modifications: catch ValueErrors and control verbosity settings

Adds option to capture ValueErrors when running batch simulation, so a single failed simulation won't halt the entire batch process. In addition, adjusts the verbose settings to capture output if verbose is false, so multiple parallel processes won't all be printing the nonlinear simulation progress bar over each other.

* Transition function edge case handling + minor stylistic edit

Handles situations in which the optimizer searches very large or small values of the 'a' parameter, which may result in a RuntimeWarning due to the exponential in the expression

* Create goodness-of-fit score files

Class file and file to hold helper functions

* Sync with laptop

* Base functional gof version

* Added gof score plotting and verbose options

* Edits to make Python gof results match Matlabs within floating pt. error

* Added batch score calculation

* Fix score array datatype, remove batch verbose printing

* Minor bug fix, imports and plotting

* Gof working update

* Added GoF example

* Added forkserver context

Added option to use the forkserver context for spawning subprocesses using the multiprocessing module.

* Updated batch simulation example 

Added forkserver example usage

* Delete examples/test-gof.ipynb

* Documentation and tests update

Updated documentation of new code to meet contribution requirements. Updated tests to reflect minor changes in code. Updated formatting to meet tox requirements.

* Update example files

Fixed filepath error, updated to match recent changes

* Documentation update

Added relevant files to docs/

* Added PyWavelets to setup

* Update build-and-deploy-docs.yml

.github/workflows/build-and-deploy-docs.yml

@@ -13,9 +13,9 @@ on:
 
 # Sets permissions of the GITHUB_TOKEN to allow deployment to GitHub Pages
 permissions:
-  contents: read
-  pages: write
   id-token: write
+  pages: write
+  contents: read
 
 # Allow only one concurrent deployment, skipping runs queued between the run in-progress and latest queued.
 # However, do NOT cancel in-progress runs as we want to allow these production deployments to complete.
@@ -38,6 +38,7 @@ jobs:
       - name: Sphinx build
         run: |
           cd docs
+          sudo apt install graphviz
           make clean html
       - name: Sanity check
         run: |

PySeismoSoil/__init__.py

@@ -1,3 +1,3 @@
 # Author: Jian Shi
 
-__version__ = 'v0.5.4'
+__version__ = 'v0.6.0'

PySeismoSoil/class_Vs_profile.py

@@ -32,6 +32,9 @@ class Vs_Profile:
          |      ...      |   ...    |   ...   |   ...   |      ...        |
          +---------------+----------+---------+---------+-----------------+
 
+        Note that the surface layer is recommended to have a thickness of at
+        least 1.0 meter, to ensure realistic overburden pressure.
+
         (The "material numbers" are integer indices that map each layer to
         their G/Gmax and damping curves.)
 
@@ -121,6 +124,14 @@ def __init__(
             else:
                 material_number = np.arange(1, n_layer_tmp + 1)
 
+            # If surface layer has a thickness less than 1.0 meter
+            if thk[0] < 1.0:
+                print(
+                    'Warning in initializing Vs_Profile: surface layer '
+                    f'thickness lower than 1.0 m (user provided = {thk[0]}).',
+                    'May result in unrealistic surface layer overburden pressure.',
+                )
+
             full_data = np.column_stack((thk, vs, xi, rho, material_number))
         elif n_col == 5:
             xi = data_[:, 2]
@@ -161,7 +172,7 @@ def __init__(
             half_space = last_row.copy()
             half_space[0] = 0  # set thickness to 0 meters
             half_space[4] = 0  # set material number to 0
-            full_data = np.row_stack((full_data, half_space))
+            full_data = np.vstack((full_data, half_space))
             thk, vs, xi, rho, material_number = full_data.T
 
         self._thk = thk

PySeismoSoil/class_batch_simulation.py

@@ -24,6 +24,12 @@ class Batch_Simulation:
         A list of simulation objects. Valid simulation objects include objects
         from these classes: ``Linear_Simulation``, ``Equiv_Linear_Simulation``,
         and ``Nonlinear_Simulation``.
+    use_ctx : bool
+        For unix systems, provides the option to use the forkserver context for spawning
+        subprocesses when running simulations in batch. The forkserver context is recommended
+        to avoid slowdowns when PySeismoSoil is being run in batch as part of a code that
+        contains additional non-PySeismoSoil variables and module imports. If use_ctx is
+        set to True, the top-level code must be guarded under `if __name__ == "__main__":`.
 
     Attributes
     ----------
@@ -42,7 +48,11 @@ class Batch_Simulation:
         When ``list_of_simulations`` has length 0
     """
 
-    def __init__(self, list_of_simulations: list[Simulation_Results]) -> None:
+    def __init__(
+            self,
+            list_of_simulations: list[Simulation_Results],
+            use_ctx: bool = False,
+    ) -> None:
         if not isinstance(list_of_simulations, list):
             raise TypeError('`list_of_simulations` should be a list.')
 
@@ -73,11 +83,23 @@ def __init__(self, list_of_simulations: list[Simulation_Results]) -> None:
         self.n_simulations = n_simulations
         self.sim_type = type(sim_0)
 
+        self.use_ctx = use_ctx
+        if use_ctx:
+            self.ctx = mp.get_context('forkserver')
+            self.ctx.set_forkserver_preload([
+                'PySeismoSoil.class_Vs_profile',
+                'PySeismoSoil.class_ground_motion',
+                'PySeismoSoil.class_simulation',
+                'PySeismoSoil.class_batch_simulation',
+            ])
+
     def run(
             self,
             parallel: bool = False,
             n_cores: int | None = 1,
             base_output_dir: str | None = None,
+            catch_errors: bool = False,
+            verbose: bool = True,
             options: dict[str, Any] | None = None,
     ) -> list[Simulation_Results]:
         """
@@ -93,6 +115,13 @@ def run(
         base_output_dir : str | None
             The parent directory for saving the output files/figures of the
             current batch.
+        catch_errors : bool
+            Optionally allows for ValueErrors to be caught during batch simulation,
+            so a single error simulation doesn't interrupt the running of others in the
+            batch. Simulations that have caught errors will be replaced by `None` in the
+            results list.
+        verbose : bool
+            Whether to print the parallel computing progress info.
         options : dict[str, Any] | None
             Options to be passed to the ``run()`` methods of the relevant
             simulation classes (linear, equivalent linear, or nonlinear). Check
@@ -114,24 +143,40 @@ def run(
             current_time = hlp.get_current_time(for_filename=True)
             base_output_dir = os.path.join('./', 'batch_sim_%s' % current_time)
 
-        other_params = [n_digits, base_output_dir, options]
+        other_params = [n_digits, base_output_dir, catch_errors, options]
 
         if not parallel:
             sim_results = []
             for i in range(self.n_simulations):
                 sim_results.append(self._run_single_sim([i, other_params]))
             # END FOR
         else:
-            # Because no outputs can be printed to stdout in the parellel pool
-            if options.get('verbose', True):  # default value is `True`
-                print('Parallel computing in progress...', end=' ')
-
-            p = mp.Pool(n_cores)
-            sim_results = p.map(
-                self._run_single_sim,
-                itertools.product(range(N), [other_params]),
-            )
-            if options.get('verbose', True):
+            sim_results = []
+
+            if self.use_ctx:
+                if verbose:
+                    print(
+                        'Parallel computing in progress using forkserver...',
+                        end=' ',
+                    )
+
+                with self.ctx.Pool(processes=n_cores) as p:
+                    sim_results = p.map(
+                        self._run_single_sim,
+                        itertools.product(range(N), [other_params]),
+                    )
+
+            else:
+                if verbose:
+                    print('Parallel computing in progress...', end=' ')
+
+                with mp.Pool(n_cores) as p:
+                    sim_results = p.map(
+                        self._run_single_sim,
+                        itertools.product(range(N), [other_params]),
+                    )
+
+            if verbose:
                 print('done.')
 
             # Because no figures can be plotted in the parallel pool:
@@ -153,13 +198,14 @@ def _run_single_sim(self, all_params: list[Any]) -> Simulation_Results:
         all_params : list[Any]
             All the parameters needed for running the simulation. It should
             have the following structure:
-                [i, [n_digits, base_output_dir, options]]
+                [i, [n_digits, base_output_dir, catch_errors, options]]
             where:
                 - ``i`` is the index of the current simulation in the batch.
                 - ``n_digits`` is the number of digits of the length of the
                   batch. (For example, if there are 125 simulations, then
                   ``n_digits`` should be 3.)
                 - ``base_output_dir``: same as in the ``run()`` method
+                - ``catch_errors``: same as in the ``run()`` method
                 - ``options``: same as in the ``run()`` method
 
         Returns
@@ -168,13 +214,21 @@ def _run_single_sim(self, all_params: list[Any]) -> Simulation_Results:
             Simulation results of a single simulation object.
         """
         i, other_params = all_params  # unpack
-        n_digits, base_output_dir, options = other_params  # unpack
+        n_digits, base_output_dir, catch_errors, options = other_params  # unpack
         output_dir = os.path.join(base_output_dir, str(i).rjust(n_digits, '0'))
         if self.sim_type == Nonlinear_Simulation:
             options.update({'sim_dir': output_dir})
         else:  # linear or equivalent linear
             options.update({'output_dir': output_dir})
 
         sim_obj = self.list_of_simulations[i]
-        sim_result = sim_obj.run(**options)
+        if catch_errors:
+            try:
+                sim_result = sim_obj.run(**options)
+            except ValueError:
+                sim_result = None
+                print('Warning: ValueError encountered.')
+        else:
+            sim_result = sim_obj.run(**options)
+
         return sim_result

PySeismoSoil/class_curves.py

@@ -3,7 +3,6 @@
 import os
 from typing import Any, Literal, Type
 
-import matplotlib.pyplot as plt
 import numpy as np
 from matplotlib.axes import Axes
 from matplotlib.figure import Figure
@@ -719,9 +718,11 @@ def plot(
         plot_interpolated : bool
             Whether to plot the interpolated curve or the raw data
         fig : Figure | None
-            Figure object. If None, a new figure will be created.
+            Figure object. If None, one new figure will be created.
         ax : Axes | None
-            Axes object. If None, a new axes will be created.
+            Axes object. If None, one new Axes object will be created.
+            If the user provides an Axes object but no Figure object, the
+            user provided Axes object will be overwritten by a new Axes object.
         title : str | None
             Title of plot.
         xlabel : str | None
@@ -744,8 +745,15 @@ def plot(
         ax : Axes
             The axes object being created or being passed into this function.
         """
-        fig = plt.figure()
-        ax = plt.axes()
+        if fig is None:  # User provided ax but not fig, or user provided neither
+            fig, ax = hlp._process_fig_ax_objects(
+                fig, None, figsize=figsize, dpi=dpi
+            )
+        elif fig is not None and ax is None:  # User provided fig but not ax
+            fig, ax = hlp._process_fig_ax_objects(
+                fig, ax, figsize=figsize, dpi=dpi
+            )
+
         for curve in self.curves:
             curve.plot(
                 plot_interpolated=plot_interpolated,