FEAT: Extend extreme method to accept absolute value extremes (#51)

Author: vegard-solum-4ss

docs/user_guide/concepts_utils/directional_spectrum.rst

@@ -80,12 +80,18 @@ Calculate the mean zero-crossing period, Tz:
     spectrum.tz
 
 Calculate extreme values using the :meth:`~waveresponse.DirectionalSpectrum.extreme`
-method. The method takes two arguments: the duration of the process (in seconds),
-and a quantile, ``q``.
+method. The method takes three arguments: the duration of the process (in seconds),
+the quantile, ``q``, and a boolean flag, ``absmax``, determining whether to compute absolute
+value extremes (or only consider the maxima (`default`)).
 
 .. code-block:: python
 
     duration = 3 * 3600   # 3 hours
 
+    # Extreme maximum
     mpm = spectrum.extreme(duration, q=0.37)   # most probable maximum (MPM)
     q90 = spectrum.extreme(duration, q=0.90)   # 90-th quantile
+
+    # Extreme absolute value maximum (i.e., minima are taken into account)
+    mpm = spectrum.extreme(duration, q=0.37, absmax=True)   # most probable maximum (MPM)
+    q90 = spectrum.extreme(duration, q=0.90, absmax=True)   # 90-th quantile

docs/user_guide/standardized_spectra.rst

@@ -173,7 +173,7 @@ Torsethaugen
 ------------
 The *Torsethaugen* spectrum allows you to set up a double-peaked spectrum that represents
 sea states which includes both a remotely generated swell component (with low frequency energy)
-and a locally generated wind component (with high frequency energy). The spectral spectral model
+and a locally generated wind component (with high frequency energy). The spectral model
 was developed based on average measured spectra for Norwegian waters (Haltenbanken and Statfjord).
 The Torsethaugen spectrum is described by two parameter (i.e. :math:`H_s` and :math:`T_p`), and
 is given by:

src/waveresponse/_core.py

@@ -1471,7 +1471,7 @@ def tz(self):
         m2 = self.moment(2, freq_hz=True)
         return np.sqrt(m0 / m2)
 
-    def extreme(self, t, q=0.37):
+    def extreme(self, t, q=0.37, absmax=False):
         """
         Compute the q-th quantile extreme value (assuming a Gaussian process).
 
@@ -1483,24 +1483,28 @@ def extreme(self, t, q=0.37):
         of the process, and ``q`` is the quantile. Setting ``q=0.37`` yields the
         most probable maximum (MPM).
 
-        The extreme value indicates a level which the maximum value of the process
-        amplitudes will be below with a certain probability. Note that the method
-        only computes extreme values for the maxima, and does not consider the minima.
-
         Parameters
         ----------
         t : float
             Time/duration in seconds for which the of the process is observed.
         q : float or array-like
             Quantile or sequence of quantiles to compute. Must be between 0 and 1
             (inclusive).
+        absmax : bool
+            Whether to compute absolute value extremes (i.e., taking the minima into account).
+            If ``False`` (default), only the maxima are considered. See Notes.
 
         Returns
         -------
         x : float or array
-            Extreme value(s). During the given time period, the maximum value of
-            the process amplitudes will be below the returned value with a given
-            probability.
+            Extreme value(s). During the given time period, the maximum value (or
+            absolute value maximum) of the process amplitudes will be below the
+            returned value with the given probability.
+
+        Notes
+        -----
+        Computing absolute value extremes by setting ``absmax=True`` is equivalent
+        to doubling the expected zero-crossing rate, ``fz = 1 / Tz``.
 
         Notes
         -----
@@ -1512,8 +1516,14 @@ def extreme(self, t, q=0.37):
            Cambridge University Press.
 
         """
+
+        if absmax:
+            tz = self.tz / 2.0
+        else:
+            tz = self.tz
+
         q = np.asarray_chkfinite(q)
-        return self.std() * np.sqrt(2.0 * np.log((t / self.tz) / np.log(1.0 / q)))
+        return self.std() * np.sqrt(2.0 * np.log((t / tz) / np.log(1.0 / q)))
 
 
 class WaveSpectrum(DirectionalSpectrum):

tests/test_core.py

@@ -3374,13 +3374,13 @@ def test_extreme_float(self):
         vals = np.ones((len(freq), len(dirs)))
         spectrum = wr.DirectionalSpectrum(freq, dirs, vals, freq_hz=True, degrees=True)
 
-        sigma = spectrum.std()
-        tz = spectrum.tz
-
         T = 360 * 24 * 60.0**2
         q = 0.99
         extreme_out = spectrum.extreme(T, q=q)
 
+        sigma = spectrum.std()
+        tz = spectrum.tz
+
         extreme_expect = sigma * np.sqrt(2.0 * np.log((T / tz) / np.log(1.0 / q)))
 
         assert extreme_out == pytest.approx(extreme_expect)
@@ -3394,13 +3394,13 @@ def test_extreme_list(self):
         vals = np.ones((len(freq), len(dirs)))
         spectrum = wr.DirectionalSpectrum(freq, dirs, vals, freq_hz=True, degrees=True)
 
-        sigma = spectrum.std()
-        tz = spectrum.tz
-
         T = 360 * 24 * 60.0**2
         q = [0.1, 0.5, 0.99]
         extreme_out = spectrum.extreme(T, q=q)
 
+        sigma = spectrum.std()
+        tz = spectrum.tz
+
         extreme_expect = [
             sigma * np.sqrt(2.0 * np.log((T / tz) / np.log(1.0 / q[0]))),
             sigma * np.sqrt(2.0 * np.log((T / tz) / np.log(1.0 / q[1]))),
@@ -3418,16 +3418,38 @@ def test_extreme_mpm(self):
         vals = np.ones((len(freq), len(dirs)))
         spectrum = wr.DirectionalSpectrum(freq, dirs, vals, freq_hz=True, degrees=True)
 
-        sigma = spectrum.std()
-        tz = spectrum.tz
-
         T = 360 * 24 * 60.0**2
         extreme_out = spectrum.extreme(T, q=0.37)
 
+        sigma = spectrum.std()
+        tz = spectrum.tz
+
         extreme_expect = sigma * np.sqrt(2.0 * np.log(T / tz))
 
         assert extreme_out == pytest.approx(extreme_expect, rel=1e-3)
 
+    def test_extreme_absmax(self):
+        f0 = 0.0
+        f1 = 2.0
+
+        freq = np.linspace(f0, f1, 20)
+        dirs = np.arange(5, 360, 10)
+        vals = np.ones((len(freq), len(dirs)))
+        spectrum = wr.DirectionalSpectrum(freq, dirs, vals, freq_hz=True, degrees=True)
+
+        T = 360 * 24 * 60.0**2
+        q = 0.99
+        extreme_out = spectrum.extreme(T, q=q, absmax=True)
+
+        sigma = spectrum.std()
+        tz_absmax = spectrum.tz / 2.0
+
+        extreme_expect = sigma * np.sqrt(
+            2.0 * np.log((T / tz_absmax) / np.log(1.0 / q))
+        )
+
+        assert extreme_out == pytest.approx(extreme_expect)
+
 
 class Test_WaveSpectrum:
     def test__init__(self):