diff --git a/.github/workflows/wheelbuild-benchmark-test.yml b/.github/workflows/wheelbuild-benchmark-test.yml new file mode 100644 index 00000000..431284a4 --- /dev/null +++ b/.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 diff --git a/.gitignore b/.gitignore index 9ddc022b..481f7409 100644 --- a/.gitignore +++ b/.gitignore @@ -61,4 +61,5 @@ coverage_report* .hypothesis .mutmut-cache -.python-version \ No newline at end of file +.python-version +/load_signal.csv diff --git a/benchmarks/conftest.py b/benchmarks/conftest.py new file mode 100644 index 00000000..7dcd2ea7 --- /dev/null +++ b/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 diff --git a/benchmarks/generate_time_signal.py b/benchmarks/generate_time_signal.py new file mode 100644 index 00000000..1986ce28 --- /dev/null +++ b/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) diff --git a/benchmarks/test_rainflow_counters.py b/benchmarks/test_rainflow_counters.py new file mode 100644 index 00000000..2f0b4759 --- /dev/null +++ b/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." diff --git a/pyproject.toml b/pyproject.toml index 2c63dbb2..790121c6 100644 --- a/pyproject.toml +++ b/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] diff --git a/setup.py b/setup.py index 0856b20d..d1ceb12a 100644 --- a/setup.py +++ b/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, " diff --git a/src/pylife/stress/rainflow/extension.pyx b/src/pylife/stress/rainflow/extension.pyx new file mode 100644 index 00000000..19d4e3b4 --- /dev/null +++ b/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] diff --git a/src/pylife/stress/rainflow/fourpoint.py b/src/pylife/stress/rainflow/fourpoint.py index 5f2d2ea4..628cf5b5 100644 --- a/src/pylife/stress/rainflow/fourpoint.py +++ b/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) diff --git a/src/pylife/stress/rainflow/general.py b/src/pylife/stress/rainflow/general.py index d9b335c3..3395a5c0 100644 --- a/src/pylife/stress/rainflow/general.py +++ b/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 diff --git a/src/pylife/stress/rainflow/recorders.py b/src/pylife/stress/rainflow/recorders.py index e560d83f..0496fc32 100644 --- a/src/pylife/stress/rainflow/recorders.py +++ b/src/pylife/stress/rainflow/recorders.py @@ -30,8 +30,8 @@ class LoopValueRecorder(AbstractRecorder): def __init__(self): """Instantiate a LoopRecorder.""" super().__init__() - self._values_from = [] - self._values_to = [] + self._values_from = np.zeros((0,)) + self._values_to = np.zeros((0,)) @property def values_from(self): @@ -53,8 +53,8 @@ def collective(self): def record_values(self, values_from, values_to): """Record the loop values.""" - self._values_from += values_from - self._values_to += values_to + self._values_from = np.append(self._values_from, values_from) + self._values_to = np.append(self._values_to, values_to) def histogram_numpy(self, bins=10): """Calculate a histogram of the recorded values into a plain numpy.histogram2d. @@ -109,8 +109,8 @@ class FullRecorder(LoopValueRecorder): def __init__(self): """Instantiate a FullRecorder.""" super().__init__() - self._index_from = [] - self._index_to = [] + self._index_from = np.array([], dtype=np.uintp) + self._index_to = np.array([], dtype=np.uintp) @property def index_from(self): @@ -137,5 +137,9 @@ def collective(self): def record_index(self, index_from, index_to): """Record the index.""" - self._index_from += index_from - self._index_to += index_to + self._index_from = np.concatenate( + (self._index_from, np.asarray(index_from, dtype=np.uintp)) + ) + self._index_to = np.concatenate( + (self._index_to, np.asarray(index_to, dtype=np.uintp)) + ) diff --git a/src/pylife/stress/rainflow/threepoint.py b/src/pylife/stress/rainflow/threepoint.py index bf6e07f0..010b6b59 100644 --- a/src/pylife/stress/rainflow/threepoint.py +++ b/src/pylife/stress/rainflow/threepoint.py @@ -17,8 +17,8 @@ __author__ = "Johannes Mueller" __maintainer__ = __author__ -import cython import numpy as np +from pylife.rainflow_ext import threepoint_loop from .general import AbstractDetector @@ -102,11 +102,6 @@ def __init__(self, recorder): """ super().__init__(recorder) - @cython.locals( - start=cython.int, front=cython.int, back=cython.int, - highest_front=cython.int, lowest_front=cython.int, - start_val=cython.double, front_val=cython.double, back_val=cython.double, - turns=cython.double[:]) def process(self, samples): """Process a sample chunk. @@ -124,54 +119,29 @@ def process(self, samples): if len(self._residuals) == 0: residuals = samples[:1] - residual_index = [0, 1] else: residuals = self._residuals[:-1] - residual_index = [*range(len(residuals))] 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 + turns = np.concatenate((residuals, turns_values, samples[-1:])) + turns_index = np.concatenate((self._residual_index, turns_index.astype(np.uintp))) highest_front = np.argmax(residuals) lowest_front = np.argmin(residuals) - from_vals = [] - to_vals = [] - from_index = [] - to_index = [] - - back = residual_index[-1] + 1 - while back < turns.shape[0]: - if len(residual_index) >= 2: - start = residual_index[-2] - front = residual_index[-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 - np.abs(back_val - front_val) >= np.abs(front_val - start_val)): - from_vals.append(start_val) - to_vals.append(front_val) - from_index.append(turns_index[start]) - to_index.append(turns_index[front]) - residual_index.pop() - residual_index.pop() - continue - - residual_index.append(back) - back += 1 + ( + from_vals, + to_vals, + from_index, + to_index, + residual_index + ) = threepoint_loop(turns, turns_index, highest_front, lowest_front, len(residuals)) self._recorder.record_values(from_vals, to_vals) self._recorder.record_index(from_index, to_index) - self._residuals = turns_np[residual_index] + self._residuals = turns[residual_index] self._residual_index = turns_index[residual_index[:-1]] self._recorder.report_chunk(len(samples)) diff --git a/tests/stress/rainflow/test_recorders.py b/tests/stress/rainflow/test_recorders.py index a765a558..7393e97c 100644 --- a/tests/stress/rainflow/test_recorders.py +++ b/tests/stress/rainflow/test_recorders.py @@ -123,8 +123,8 @@ def test_full_rainflow_recorder_empty_collective_default(): expected = pd.DataFrame({ 'from': [], 'to': [], - 'index_from': [], - 'index_to': [] + 'index_from': pd.Series([], dtype=np.uintp), + 'index_to': pd.Series([], dtype=np.uintp) }) pd.testing.assert_frame_equal(fr.collective, expected) @@ -140,10 +140,12 @@ def test_full_rainflow_recorder_two_non_zero_collective(): expected = pd.DataFrame({ 'from': [vf1, vf2], 'to': [vt1, vt2], - 'index_from': [if1, if2], - 'index_to': [it1, it2] + 'index_from': pd.Series([if1, if2], dtype=np.uintp), + 'index_to': pd.Series([it1, it2], dtype=np.uintp) }) + print(expected) + print(expected.dtypes) pd.testing.assert_frame_equal(fr.collective, expected) diff --git a/tests/stress/rainflow/test_threepoint.py b/tests/stress/rainflow/test_threepoint.py index 9739254b..55b26fa6 100644 --- a/tests/stress/rainflow/test_threepoint.py +++ b/tests/stress/rainflow/test_threepoint.py @@ -381,6 +381,44 @@ def test_residual_index(self): np.testing.assert_array_equal(self._dtor.residual_index, np.array([0, 1, 6])) +class TestThreePointDampeningSplit(unittest.TestCase): + r''' + 1 2 3 4 5 6 + ------------------------------------------------------------------------------------------------------- + 6 R | | | | | | | + -----------/-\----------------------------------------------------------------------------------------- + 5 / \ 2 | | | | | | | + ---------/-----\---------/-\--------------------------------------------------------------------------- + 4 / \ / \ 1 | | | 1 | | | | + -------/---------\-----/-----\-/-\--------------------------------------------------------------------- + 3 / \ / 1---\ | | 1 | | | | | + -----/-------------\-/-------------\------------------------------------------------------------------- + 2 / 2---------------\ | | | | | | | + ---/---------------------------------\----------------------------------------------------------------- + 1 R R | | | | | | | + ------------------------------------------------------------------------------------------------------- + 0 1 2 3 4 5 6 + ''' + def setUp(self): + signal = np.array([1., 6., 2., 5., 3., 4., 1.]) + self._fr, self._dtor = process_signal(signal[:5]) + self._dtor.process(signal[5:]) + + def test_values(self): + np.testing.assert_array_equal(self._fr.values_from, np.array([3., 2.])) + np.testing.assert_array_equal(self._fr.values_to, np.array([4., 5.])) + + def test_indeces(self): + np.testing.assert_array_equal(self._fr.index_from, np.array([4, 2])) + np.testing.assert_array_equal(self._fr.index_to, np.array([5, 3])) + + def test_residuals(self): + np.testing.assert_array_equal(self._dtor.residuals, np.array([1., 6., 1.])) + + def test_residual_index(self): + np.testing.assert_array_equal(self._dtor.residual_index, np.array([0, 1, 6])) + + test_signals = [ np.array([0., 1., 0., 1., -1., 1., 0., 1., -1., 1., 0]), np.array([2., 5., 3., 6., 2., 3., 1., 6., 1., 4., 2., 3., 1., 4., 2., 5., 3., 4., 2.]), diff --git a/tools/odbclient/pyproject.toml b/tools/odbclient/pyproject.toml index 98ce6044..9f53a7de 100644 --- a/tools/odbclient/pyproject.toml +++ b/tools/odbclient/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"] build-backend = "setuptools.build_meta" [tool.setuptools_scm]