Introduce upcrossing analysis module (#252)

* Add test for loads.extreme.global_peaks function

The loads_extreme.global_peaks function was previously missing a test.

The test uses a simple function which can be independently analysed. The results of global_peaks and the independent analysis are then compared.

* Introduce upcrossing module

Previously there was no general means of performing an upcrossing analysis. The load.extreme.global_peaks function could only calculate peaks.

The module provides some common methods but also the ability for the user to define their own function over the zero crossing points.

* Implement loads.extreme.global_peaks in terms of upcrossing module

With the recent addition of the upcrossing module, we can implement the loads.extreme.global_peaks function using it.

* minor linting, module docstring, update parameter validation

* fix upcrossing docstring

* move upcrossing module to utils

* update upcrossing import in example

* fix last upcrossing docstring and move test to test/utils

* update description in the upcrossing notebook

* update import of upcrossing into test_upcrossing

* typo in test file

* final typo fix

---------

Co-authored-by: akeeste <[email protected]>

mhkit/loads/extreme.py

@@ -2,6 +2,7 @@
 import pandas as pd
 from scipy import stats, optimize, signal
 from mhkit.wave.resource import frequency_moment
+from mhkit.utils import upcrossing, custom
 
 def _peaks_over_threshold(peaks, threshold, sampling_rate):
     threshold_unit = np.percentile(peaks, 100*threshold, method='hazen')
@@ -61,24 +62,23 @@ def global_peaks(t, data):
     assert isinstance(t, np.ndarray), 't must be of type np.ndarray'
     assert isinstance(data, np.ndarray), 'data must be of type np.ndarray'
 
-    # eliminate zeros
-    zeroMask = (data == 0)
-    data[zeroMask] = 0.5 * np.min(np.abs(data))
-    # zero up-crossings
-    diff = np.diff(np.sign(data))
-    zeroUpCrossings_mask = (diff == 2) | (diff == 1)
-    zeroUpCrossings_index = np.where(zeroUpCrossings_mask)[0]
-    zeroUpCrossings_index = np.append(zeroUpCrossings_index, len(data) - 1)
-    # global peaks
-    npeaks = len(zeroUpCrossings_index)
-    peaks = np.array([])
-    t_peaks = np.array([])
-    for i in range(npeaks - 1):
-        peak_index = np.argmax(
-            data[zeroUpCrossings_index[i]:zeroUpCrossings_index[i + 1]])
-        t_peaks = np.append(t_peaks, t[zeroUpCrossings_index[i] + peak_index])
-        peaks = np.append(peaks, data[zeroUpCrossings_index[i] + peak_index])
-    return t_peaks, peaks
+    # Find zero up-crossings
+    inds = upcrossing(t, data)
+
+    # We also include the final point in the dataset
+    inds = np.append(inds, len(data)-1)
+
+    # As we want to return both the time and peak
+    # values, look for the index at the peak.
+    # The call to argmax gives us the index within the
+    # upcrossing period. Therefore to get the index in the
+    # original array we need to add on the index that
+    # starts the zero crossing period, ind1.
+    func = lambda ind1, ind2: np.argmax(data[ind1:ind2]) + ind1
+
+    peak_inds = np.array(custom(t, data, func, inds), dtype=int)
+
+    return t[peak_inds], data[peak_inds]
 
 
 def number_of_short_term_peaks(n, t, t_st):

mhkit/tests/loads/test_extreme.py

@@ -0,0 +1,51 @@
+import numpy as np
+import unittest
+import mhkit.loads as loads
+from numpy.testing import assert_allclose
+
+
+class TestExtreme(unittest.TestCase):
+    @classmethod
+    def setUpClass(self):
+        self.t, self.signal = self._example_waveform(self)
+
+    def _example_waveform(self):
+        # Create simple wave form to analyse.
+        # This has been created to perform
+        # a simple independent calcuation that
+        # the mhkit functions can be tested against.
+
+        A = np.array([0.5, 0.6, 0.3])
+        T = np.array([3, 2, 1])
+        w = 2 * np.pi / T
+
+        t = np.linspace(0, 4.5, 100)
+
+        signal = np.zeros(t.size)
+        for i in range(A.size):
+            signal += A[i] * np.sin(w[i] * t)
+
+        return t, signal
+
+    def _example_crest_analysis(self, t, signal):
+        # NB: This only works due to the construction
+        # of our test signal. It is not suitable as
+        # a general approach.
+        grad = np.diff(signal)
+
+        # +1 to get the index at turning point
+        turning_points = np.flatnonzero(grad[1:] * grad[:-1] < 0) + 1
+
+        crest_inds = turning_points[signal[turning_points] > 0]
+        crests = signal[crest_inds]
+
+        return crests, crest_inds
+
+    def test_global_peaks(self):
+        peaks_t, peaks_val = loads.extreme.global_peaks(self.t, self.signal)
+
+        test_crests, test_crests_ind = self._example_crest_analysis(
+            self.t, self.signal)
+
+        assert_allclose(peaks_t, self.t[test_crests_ind])
+        assert_allclose(peaks_val, test_crests)

mhkit/tests/utils/test_upcrossing.py

@@ -0,0 +1,142 @@
+from mhkit.utils import upcrossing, peaks, troughs, heights, periods, custom
+import unittest
+from numpy.testing import assert_allclose
+import numpy as np
+from scipy.optimize import fsolve
+
+
+class TestUpcrossing(unittest.TestCase):
+    @classmethod
+    def setUpClass(self):
+        self.t = np.linspace(0, 4, 1000)
+
+        self.signal = self._example_waveform(self, self.t)
+
+        # Approximiate points for the zero crossing,
+        # used as starting points in numerical
+        # solution.
+        self.zero_cross_approx = [0, 2.1, 3, 3.8]
+
+    def _example_waveform(self, t):
+        # Create simple wave form to analyse.
+        # This has been created to perform
+        # a simple independent calcuation that
+        # the mhkit functions can be tested against.
+
+        A = np.array([0.5, 0.6, 0.3])
+        T = np.array([3, 2, 1])
+        w = 2 * np.pi / T
+
+        signal = np.zeros(t.size)
+        for i in range(A.size):
+            signal += A[i] * np.sin(w[i] * t)
+
+        return signal
+
+    def _example_analysis(self, t, signal):
+        # NB: This only works due to the construction
+        # of our test signal. It is not suitable as
+        # a general approach.
+        grad = np.diff(signal)
+
+        # +1 to get the index at turning point
+        turning_points = np.flatnonzero(grad[1:] * grad[:-1] < 0) + 1
+
+        crest_inds = turning_points[signal[turning_points] > 0]
+        trough_inds = turning_points[signal[turning_points] < 0]
+
+        crests = signal[crest_inds]
+        troughs = signal[trough_inds]
+
+        heights = crests - troughs
+
+        zero_cross = fsolve(self._example_waveform, self.zero_cross_approx)
+        periods = np.diff(zero_cross)
+
+        return crests, troughs, heights, periods
+
+    def test_peaks(self):
+        want, _, _, _ = self._example_analysis(self.t, self.signal)
+
+        got = peaks(self.t, self.signal)
+
+        assert_allclose(got, want)
+
+    def test_troughs(self):
+        _, want, _, _ = self._example_analysis(self.t, self.signal)
+
+        got = troughs(self.t, self.signal)
+
+        assert_allclose(got, want)
+
+    def test_heights(self):
+        _, _, want, _ = self._example_analysis(self.t, self.signal)
+
+        got = heights(self.t, self.signal)
+
+        assert_allclose(got, want)
+
+    def test_periods(self):
+        _, _, _, want = self._example_analysis(self.t, self.signal)
+
+        got = periods(self.t, self.signal)
+
+        assert_allclose(got, want, rtol=1e-3, atol=1e-3)
+
+    def test_custom(self):
+        want, _, _, _ = self._example_analysis(self.t, self.signal)
+
+        # create a similar function to finding the peaks
+        def f(ind1, ind2): return np.max(self.signal[ind1:ind2])
+
+        got = custom(self.t, self.signal, f)
+
+        assert_allclose(got, want)
+
+    def test_peaks_with_inds(self):
+        want, _, _, _ = self._example_analysis(self.t, self.signal)
+
+        inds = upcrossing(self.t, self.signal)
+
+        got = peaks(self.t, self.signal, inds)
+
+        assert_allclose(got, want)
+
+    def test_trough_with_inds(self):
+        _, want, _, _ = self._example_analysis(self.t, self.signal)
+
+        inds = upcrossing(self.t, self.signal)
+
+        got = troughs(self.t, self.signal, inds)
+
+        assert_allclose(got, want)
+
+    def test_heights_with_inds(self):
+        _, _, want, _ = self._example_analysis(self.t, self.signal)
+
+        inds = upcrossing(self.t, self.signal)
+
+        got = heights(self.t, self.signal, inds)
+
+        assert_allclose(got, want)
+
+    def test_periods_with_inds(self):
+        _, _, _, want = self._example_analysis(self.t, self.signal)
+
+        inds = upcrossing(self.t, self.signal)
+
+        got = periods(self.t, self.signal, inds)
+
+        assert_allclose(got, want, rtol=1e-3, atol=1e-3)
+
+    def test_custom_with_inds(self):
+        want, _, _, _ = self._example_analysis(self.t, self.signal)
+
+        inds = upcrossing(self.t, self.signal)
+
+        # create a similar function to finding the peaks
+        def f(ind1, ind2): return np.max(self.signal[ind1:ind2])
+
+        got = custom(self.t, self.signal, f, inds)
+
+        assert_allclose(got, want)

mhkit/utils/__init__.py

@@ -1,5 +1,6 @@
 from .time_utils import matlab_to_datetime, excel_to_datetime, index_to_datetime
 from .stat_utils import get_statistics, vector_statistics, unwrap_vector, magnitude_phase, unorm
 from .cache import handle_caching, clear_cache
+from .upcrossing import upcrossing, peaks, troughs, heights, periods, custom
 
 _matlab = False  # Private variable indicating if mhkit is run through matlab