Merge pull request #99 from boschresearch/38-cythonize-rainflow-recor…

…ders

.github/workflows/wheelbuild-benchmark-test.yml

@@ -0,0 +1,62 @@
+name: Build wheels and perform benchmarks
+
+on:
+  pull_request:
+    branches: [develop, master]
+
+env:
+  CIBW_BUILD: cp37-* cp38-* cp39-* cp310-* cp311-* cp312-*
+  CIBW_MANYLINUX_X86_64_IMAGE: manylinux2014
+  CIBW_MANYLINUX_I686_IMAGE: manylinux2014
+
+jobs:
+  build_wheels:
+    name: Build wheels on for various systems
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest] #, windows-latest]
+
+    steps:
+      - uses: actions/checkout@v3
+        with:
+          fetch-depth: 0
+
+      - uses: actions/setup-python@v4
+        name: Install Python
+        with:
+          python-version: '3.12'
+
+      - name: Build wheels
+        run: |
+          pip wheel --no-deps -w dist .
+      - uses: actions/upload-artifact@v4
+        with:
+          name: wheel
+          path: dist/
+
+  benchmark:
+    name: Benchmark tests
+    needs: [build_wheels]
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v3
+        with:
+          fetch-depth: 0
+      - uses: actions/download-artifact@v4
+        with:
+          name: wheel
+          path: dist/
+      - uses: actions/setup-python@v4
+        name: Install Python
+        with:
+          python-version: '3.12'
+      - name: Install wheel
+        run: pip install $(ls -1 dist/*.whl)
+      - name: Install pytest
+        run: pip install pytest pytest-cov
+      - name: Generate test signal
+        run: python benchmarks/generate_time_signal.py
+      - name: Run benchmarks
+        run: pytest --no-cov -rP benchmarks

.gitignore

@@ -61,4 +61,5 @@ coverage_report*
 
 .hypothesis
 .mutmut-cache
-.python-version
+.python-version
+/load_signal.csv

benchmarks/conftest.py

@@ -0,0 +1,25 @@
+import pytest
+
+import time
+import numpy as np
+
+
+@pytest.fixture
+def benchmarker():
+
+    def the_benchmarker(rainflow_counter):
+        load_signal = np.loadtxt('load_signal.csv')
+
+        tic = time.perf_counter()
+
+        rainflow_counter.process(load_signal)
+
+        toc = time.perf_counter()
+        elapsed = toc - tic
+
+        classname = rainflow_counter.__class__.__name__
+        print(f"Processing {classname} took {elapsed:0.4f} seconds")
+
+        return elapsed
+
+    return the_benchmarker

benchmarks/generate_time_signal.py

@@ -0,0 +1,22 @@
+
+import numpy as np
+import pylife.stress.timesignal as TS
+
+np.random.seed(23424711)
+
+load_signal = TS.TimeSignalGenerator(
+    10,
+    {
+        'number': 50_000,
+        'amplitude_median': 50.0,
+        'amplitude_std_dev': 0.5,
+        'frequency_median': 4,
+        'frequency_std_dev': 3,
+        'offset_median': 0,
+        'offset_std_dev': 0.4,
+    },
+    None,
+    None,
+).query(1_000_000)
+
+np.savetxt('load_signal.csv', load_signal)

benchmarks/test_rainflow_counters.py

@@ -0,0 +1,19 @@
+import pylife.stress.rainflow as RF
+
+
+def test_threepoint(benchmarker):
+    benchmark_time = 0.05  # seconds
+    elapsed = benchmarker(RF.ThreePointDetector(recorder=RF.FullRecorder()))
+
+    assert (
+        elapsed < benchmark_time
+    ), f"Benchmark time of {benchmark_time} s not exceeded. Needed {elapsed:0.4f} s."
+
+
+def test_fourpoint(benchmarker):
+    benchmark_time = 0.05  # seconds
+    elapsed = benchmarker(RF.FourPointDetector(recorder=RF.FullRecorder()))
+
+    assert (
+        elapsed < benchmark_time
+    ), f"Benchmark time of {benchmark_time} s not exceeded. Needed {elapsed:0.4f} s."

pyproject.toml

@@ -1,6 +1,6 @@
 [build-system]
 # AVOID CHANGING REQUIRES: IT WILL BE UPDATED BY PYSCAFFOLD!
-requires = ["setuptools>=46.1.0", "setuptools_scm[toml]>=5", "wheel"]
+requires = ["setuptools>=46.1.0", "setuptools_scm[toml]>=5", "wheel", "cython", "numpy"]
 build-backend = "setuptools.build_meta"
 
 [tool.setuptools_scm]

setup.py

@@ -7,10 +7,22 @@
     Learn more under: https://pyscaffold.org/
 """
 from setuptools import setup
+from distutils.core import Extension
+import numpy
+
+ext = Extension(
+    name='pylife.rainflow_ext',
+    sources=['src/pylife/stress/rainflow/extension.pyx'],
+    include_dirs=[numpy.get_include()],
+    extra_compile_args=["-O3"]
+)
 
 if __name__ == "__main__":
     try:
-        setup(use_scm_version={"version_scheme": "no-guess-dev"})
+        setup(
+            use_scm_version={"version_scheme": "no-guess-dev"},
+            ext_modules=[ext]
+        )
     except:  # noqa
         print(
             "\n\nAn error occurred while building the project, "

src/pylife/stress/rainflow/extension.pyx

@@ -0,0 +1,138 @@
+
+
+cimport cython
+import numpy as np
+from libc.math cimport fabs
+
+
+@cython.boundscheck(False)  # Deactivate bounds checking
+@cython.wraparound(False)   # Deactivate negative indexing.
+def fourpoint_loop(double [::1] turns, size_t [::1] turns_index):
+    cdef Py_ssize_t len_turns = len(turns)
+
+    from_vals = np.empty(len_turns//2, dtype=np.float64)
+    to_vals = np.empty(len_turns//2, dtype=np.float64)
+    from_index = np.empty(len_turns//2, dtype=np.uintp)
+    to_index = np.empty(len_turns//2, dtype=np.uintp)
+
+    cdef double [::1] from_vals_v = from_vals
+    cdef double [::1] to_vals_v = to_vals
+    cdef size_t [::1] from_index_v = from_index
+    cdef size_t [::1] to_index_v = to_index
+
+    residual_index = np.empty(len_turns, dtype=np.uintp)
+    cdef size_t [::1] residual_index_v = residual_index
+
+    residual_index_v[0] = 0
+    residual_index_v[1] = 1
+
+    cdef size_t i = 2
+    cdef size_t ri = 2
+    cdef size_t t = 0
+
+    cdef double a
+    cdef double b
+    cdef double c
+    cdef double d
+    cdef double ab
+    cdef double bc
+    cdef double cd
+
+    while i < len_turns:
+        if ri < 3:
+            residual_index_v[ri] = i
+            ri += 1
+            i += 1
+            continue
+
+        a = turns[residual_index_v[ri-3]]
+        b = turns[residual_index_v[ri-2]]
+        c = turns[residual_index_v[ri-1]]
+        d = turns[i]
+
+        ab = fabs(a - b)
+        bc = fabs(b - c)
+        cd = fabs(c - d)
+        if bc <= ab and bc <= cd:
+            from_vals_v[t] = b
+            to_vals_v[t] = c
+
+            ri -= 1
+            to_index_v[t] = turns_index[residual_index_v[ri]]
+            ri -= 1
+            from_index_v[t] = turns_index[residual_index_v[ri]]
+            t += 1
+            continue
+
+        residual_index_v[ri] = i
+        ri += 1
+        i += 1
+
+    return from_vals[:t], to_vals[:t], from_index[:t], to_index[:t], residual_index[:ri]
+
+
+cpdef double _max(double a, double b):
+    return a if a > b else b
+
+
+@cython.boundscheck(False)  # Deactivate bounds checking
+@cython.wraparound(False)   # Deactivate negative indexing.
+def threepoint_loop(double [::1] turns, size_t [::1] turns_index, size_t highest_front,
+                    size_t lowest_front, size_t residual_length):
+    cdef Py_ssize_t len_turns = len(turns)
+
+    from_vals = np.empty(len_turns//2, dtype=np.float64)
+    to_vals = np.empty(len_turns//2, dtype=np.float64)
+    from_index = np.empty(len_turns//2, dtype=np.uintp)
+    to_index = np.empty(len_turns//2, dtype=np.uintp)
+
+    cdef double [::1] from_vals_v = from_vals
+    cdef double [::1] to_vals_v = to_vals
+    cdef size_t [::1] from_index_v = from_index
+    cdef size_t [::1] to_index_v = to_index
+
+    residual_index = np.empty(len_turns, dtype=np.uintp)
+    residual_index[:residual_length] = np.arange(residual_length)
+    cdef size_t [::1] residual_index_v = residual_index
+
+    residual_index_v[0] = 0
+    residual_index_v[1] = 1
+
+    cdef size_t ri = 2
+    cdef size_t t = 0
+
+    cdef size_t back = residual_index_v[1] + 1
+    cdef size_t front
+    cdef size_t start
+
+    cdef double start_val
+    cdef double front_val
+    cdef double back_valstar
+
+    while back < len_turns:
+        if ri >= 2:
+            start = residual_index_v[ri-2]
+            front = residual_index_v[ri-1]
+            start_val, front_val, back_val = turns[start], turns[front], turns[back]
+
+            if front_val > turns[highest_front]:
+                highest_front = front
+            elif front_val < turns[lowest_front]:
+                lowest_front = front
+            elif (start >= _max(lowest_front, highest_front) and
+                  fabs(back_val - front_val) >= fabs(front_val - start_val)):
+                from_vals_v[t] = start_val
+                to_vals_v[t] = front_val
+
+                from_index_v[t] = turns_index[start]
+                to_index_v[t] = turns_index[front]
+
+                t += 1
+                ri -= 2
+                continue
+
+        residual_index[ri] = back
+        ri += 1
+        back += 1
+
+    return from_vals[:t], to_vals[:t], from_index[:t], to_index[:t], residual_index[:ri]

src/pylife/stress/rainflow/fourpoint.py

@@ -19,6 +19,7 @@
 
 
 import numpy as np
+from pylife.rainflow_ext import fourpoint_loop
 
 from .general import AbstractDetector
 
@@ -107,43 +108,15 @@ def process(self, samples):
         turns_index, turns_values = self._new_turns(samples)
 
         turns_np = np.concatenate((residuals, turns_values, samples[-1:]))
-        turns_index = np.concatenate((self._residual_index, turns_index))
-
-        turns = turns_np
-
-        from_vals = []
-        to_vals = []
-        from_index = []
-        to_index = []
-
-        residual_index = [0, 1]
-        i = 2
-        while i < len(turns):
-            if len(residual_index) < 3:
-                residual_index.append(i)
-                i += 1
-                continue
-
-            a = turns_np[residual_index[-3]]
-            b = turns_np[residual_index[-2]]
-            c = turns_np[residual_index[-1]]
-            d = turns_np[i]
-
-            ab = np.abs(a - b)
-            bc = np.abs(b - c)
-            cd = np.abs(c - d)
-            if bc <= ab and bc <= cd:
-                from_vals.append(b)
-                to_vals.append(c)
-
-                idx_2 = turns_index[residual_index.pop()]
-                idx_1 = turns_index[residual_index.pop()]
-                from_index.append(idx_1)
-                to_index.append(idx_2)
-                continue
-
-            residual_index.append(i)
-            i += 1
+        turns_index = np.concatenate((self._residual_index, turns_index.astype(np.uintp)))
+
+        (
+            from_vals,
+            to_vals,
+            from_index,
+            to_index,
+            residual_index
+        ) = fourpoint_loop(turns_np, turns_index)
 
         self._recorder.record_values(from_vals, to_vals)
         self._recorder.record_index(from_index, to_index)

src/pylife/stress/rainflow/general.py

@@ -60,7 +60,7 @@ def find_turns(samples):
 
     def clean_nans(samples):
         nans = pd.isna(samples)
-        if any(nans):
+        if nans.any():
             warnings.warn(UserWarning("At least one NaN like value has been dropped from the input signal."))
             return samples[~nans], nans
         return samples, None
@@ -138,7 +138,7 @@ def __init__(self, recorder):
         self._sample_tail = np.array([])
         self._recorder = recorder
         self._head_index = 0
-        self._residual_index = np.array([0], dtype=np.int64)
+        self._residual_index = np.array([0], dtype=np.uintp)
         self._residuals = np.array([])
 
     @property