ENH: Add the gufunc backlash_sum (Prandtl-Ishlinskii model)

Author: WarrenWeckesser

README.md

@@ -121,6 +121,7 @@ processed with the script in `ufunclab/tools/conv_template.py`.
 | --------                                        | -----------                                           |
 | [`fillnan1d`](#fillnan1d)                       | Replace `nan` using linear interpolation              |
 | [`backlash`](#backlash)                         | Backlash operator                                     |
+| [`backlash_sum`](#backlash_sum)                 | Sum backlash operators (Prandtl-Ishlinkskii)          |
 | [`hysteresis_relay`](#hysteresis_relay)         | Relay with hysteresis (Schmitt trigger)               |
 | [`sosfilter`](#sosfilter)                       | SOS (second order sections) linear filter             |
 | [`sosfilter_ic`](#sosfilter_ic)                 | SOS linear filter with initial condition              |
@@ -1727,6 +1728,43 @@ the plot
 
 ![Backlash plot](https://github.com/WarrenWeckesser/ufunclab/blob/main/examples/backlash_demo.png)
 
+#### `backlash_sum`
+
+`backlash_sum` computes the linear combination of `m` `backlash` operations.
+This is known as the Prandtl-Ishlinskii hysteresis model.  The function is a gufunc
+with signature `(n),(m),(m),(m)->(n),(m)`.  The second return value is the final value of
+each of the backlash processes.
+
+For example,
+
+```
+>>> import numpy as np
+>>> from ufunclab import backlash_sum
+
+>>> x = np.array([0.0, 0.2, 0.5, 1.1, 1.25, 1.0, 0.2, -1])
+
+Here the weights are all the same and happen to sum to 1, but that
+is not required in general.
+
+>>> w = np.array([0.25, 0.25, 0.25, 0.25])
+>>> deadband = np.array([0.2, 0.4, 0.6, 0.8])
+>>> initial = np.zeros(4)
+
+>>> y, final = backlash_sum(x, w, deadband, initial)
+>>> y
+array([ 0.    ,  0.025 ,  0.25  ,  0.85  ,  1.    ,  0.9875,  0.45  , -0.75  ])
+```
+
+-----
+
+Another example is in the script `backlash_sum_demo.py` in the `examples` directory.
+It passes two cycles of a sine wave through a Prandtl-Ishlinskii model with three
+backlash operations.  The weights are `w = [0.25, 1.0, 0.5]`, the deadband values
+are `[0.5, 1.0, 2.0]` and the initial values are all zero.
+
+![Backlash_sum plot](https://github.com/WarrenWeckesser/ufunclab/blob/main/examples/backlash_sum_demo_x_vs_t.png)
+
+![Backlash_sum plot](https://github.com/WarrenWeckesser/ufunclab/blob/main/examples/backlash_sum_demo_y_vs_x.png)
 
 #### `hysteresis_relay`
 

examples/backlash_sum_demo.py

@@ -0,0 +1,34 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from ufunclab import backlash_sum
+
+
+t = np.linspace(0, 4*np.pi, 1000)
+x = np.sin(t)
+
+w = np.array([0.25, 1.0, 1.25])
+deadband = np.array([0.5, 1.0, 2.0])
+initial = np.zeros(3)
+
+y, final = backlash_sum(x, w, deadband, initial)
+
+
+plt.figure(figsize=(6, 4))
+plt.plot(t, x, linewidth=2, label='x(t)')
+plt.plot(t, y, '--', linewidth=2, label='backlash_sum(x(t))')
+plt.xlabel('t')
+plt.title('backlash_sum')
+plt.legend(shadow=True)
+plt.grid(alpha=0.6)
+plt.savefig('backlash_sum_demo_x_vs_t.png')
+
+plt.figure(figsize=(6, 4))
+plt.plot(x, y, linewidth=2)
+plt.xlabel('x')
+plt.ylabel('y')
+plt.title('backlash_sum')
+plt.grid(alpha=0.6)
+plt.savefig('backlash_sum_demo_y_vs_x.png')
+
+# plt.show()

examples/backlash_sum_demo_x_vs_t.png

@@ -8,7 +8,7 @@ requires = [
 
 [project]
 name = 'ufunclab'
-version = '0.0.8.dev14'
+version = '0.0.8.dev15'
 description = 'NumPy ufuncs and utilities.'
 readme = 'README.md'
 requires-python = '>=3.9'

examples/backlash_sum_demo_y_vs_x.png

@@ -42,4 +42,51 @@ backlash_core(
     }
 }
 
+
+template<typename T>
+static void
+backlash_sum_core(
+        npy_intp n,                     // core dimension n
+        npy_intp m,                     // core dimension m
+        T *p_x,                         // pointer to first element of x, a strided 1-d array with shape (n,)
+        const npy_intp x_stride,        // stride (in bytes) of x
+        T *p_w,                         // pointer to first element of w, a strided 1-d array with shape (m,)
+        const npy_intp w_stride,        // stride (in bytes) of w
+        T *p_deadband,                  // pointer to first element of deadband, a strided 1-d array with shape (m,)
+        const npy_intp deadband_stride, // stride (in bytes) of deadband
+        T *p_initial,                   // pointer to initial
+        const npy_intp initial_stride,  // stride (in bytes) of initial, a strided 1-d array with shape (m,)
+        T *p_out,                       // pointer to out, a strided 1-d array with shape (n,)
+        const npy_intp out_stride,      // stride (in bytes) of out
+        T *p_final,                     // point to final, a strided 1-d array with shape (m,)
+        const npy_intp final_stride     // stride (in bytes) of final
+)
+{
+    for (npy_intp j = 0; j < m; ++j) {
+        set(p_final, final_stride, j, get(p_initial, initial_stride, j));
+    }
+    for (npy_intp k = 0; k < n; ++k) {
+        T current_x = get(p_x, x_stride, k);
+        T y_total = 0.0;
+        for (npy_intp j = 0; j < m; ++j) {
+            T current_y = get(p_final, final_stride, j);
+            T halfband = get(p_deadband, deadband_stride, j)/2;
+            T xminus = current_x - halfband;
+            if (xminus > current_y) {
+                current_y = xminus;
+            }
+            else {
+                T xplus = current_x + halfband;
+                if (xplus < current_y) {
+                    current_y = xplus;
+                }
+            }
+            set(p_final, final_stride, j, current_y);
+            T w = get(p_w, w_stride, j);
+            y_total += w*current_y;
+        }
+        set(p_out, out_stride, k, y_total);
+    }
+}
+
 #endif

pyproject.toml

@@ -29,6 +29,49 @@
 array([0. , 0.8, 0.9, 0.9, 1.3, 1.3, 1.3, 0.7])
 """
 
+BACKLASH_SUM_DOCSTRING = """\
+backlash_sum(x, w, deadband, initial, /, ...)
+
+Compute the linear combination of several backlash processes. This
+operation is also known as the Prandtl-Ishlinskii hysteresis model.
+
+Parameters
+----------
+x : array_like, length n
+    The input signal
+w : array_like, length m
+    The weights of the backlash operators.
+deadband : array_like, length m
+    The deadband values of the backlash operators.
+initial : array_like, length m
+    The initial values of the backlash operators.
+
+Returns
+-------
+out : ndarray, same shape as `x`
+    Output of the Prandtl-Ishlinskii hysteresis model.
+final : ndarray, same shape as `initial`
+    State of the individual backlash processes.
+
+Examples
+--------
+>>> import numpy as np
+>>> from ufunclab import backlash_sum
+
+>>> x = np.array([0.0, 0.2, 0.5, 1.1, 1.25, 1.0, 0.2, -1])
+
+Here the weights are all the same and happen to sum to 1, but that
+is not required in general.
+
+>>> w = np.array([0.25, 0.25, 0.25, 0.25])
+>>> deadband = np.array([0.2, 0.4, 0.6, 0.8])
+>>> initial = np.zeros(4)
+
+>>> y, final = backlash_sum(x, w, deadband, initial)
+>>> y
+array([ 0.    ,  0.025 ,  0.25  ,  0.85  ,  1.    ,  0.9875,  0.45  , -0.75  ])
+
+"""
 
 backlash_core_source = UFuncSource(
     funcname='backlash_core',
@@ -44,12 +87,25 @@
 )
 
 
+backlash_sum_core_source = UFuncSource(
+    funcname='backlash_sum_core',
+    typesignatures=['ffff->ff', 'dddd->dd', 'gggg->gg'],
+)
+
+backlash_sum_gufunc = UFunc(
+    name='backlash_sum',
+    docstring=BACKLASH_SUM_DOCSTRING,
+    header='backlash_gufunc.h',
+    signature='(n),(m),(m),(m) -> (n),(m)',
+    sources=[backlash_sum_core_source],
+)
+
 MODULE_DOCSTRING = """\
-This module defines the backlash function.
+This module defines the backlash and backlash_sum functions.
 """
 
 extmod = UFuncExtMod(
     module='_backlash',
     docstring=MODULE_DOCSTRING,
-    ufuncs=[backlash_gufunc],
+    ufuncs=[backlash_gufunc, backlash_sum_gufunc],
 )

src/backlash/backlash_gufunc.h

@@ -62,6 +62,7 @@
     'tri_area': '._tri_area',
     'tri_area_indexed': '._tri_area',
     'backlash': '._backlash',
+    'backlash_sum': '._backlash',
     'fillnan1d': '._fillnan1d',
     'linear_interp1d': '._linear_interp1d',
     'deadzone': '._deadzone',

src/backlash/define_cxx_gufunc_extmod.py

@@ -1,6 +1,6 @@
 import numpy as np
-from numpy.testing import assert_array_equal
-from ufunclab import backlash
+from numpy.testing import assert_array_equal, assert_allclose
+from ufunclab import backlash, backlash_sum
 
 
 def test_all_zeros():
@@ -50,3 +50,18 @@ def test_out_nonstandard_strides():
     y1 = backlash(x, 1, 1)
     y2 = backlash(x, 1, 1, out=out)
     assert_array_equal(y1, y2)
+
+
+def test_backlash_sum():
+    x = np.array([0.25, 0.5, 10.0, 0.125, 0.2, 0.2, 0.2, 0.2])
+    w = np.array([0.5, 3.0, 2.0])
+    deadband = np.array([0.5, 1.0, 0.75])
+    initial = np.array([1.5, 5.0, -2.0])
+    yb = []
+    for k in range(3):
+        y = backlash(x, deadband[k], initial[k])
+        yb.append(y)
+    yb = np.array(yb)
+    yy, final = backlash_sum(x, w, deadband, initial)
+    assert_allclose(yy, w@yb, rtol=1e-14)
+    assert_allclose(final, yb[:, -1], rtol=1e-14)