8.3.7. scipy on the OpenMV Cam

The scipy submodule of ulab is a small subset of CPython’s scipy, picked for usefulness on a microcontroller. It is imported as:

from ulab import scipy as sp

The submodules are:

  • sp.signal – signal processing (one function: sosfilt).

  • sp.linalg – triangular solvers and Cholesky-factorised solver.

  • sp.optimize – 1-D root finding and minimisation.

  • sp.integrate – numerical integration (quadrature).

  • sp.special – a handful of special functions (erf, gamma, …).

If you call into a function that is not compiled into your firmware, you will get an AttributeError. dir(sp) and dir(sp.signal) etc. tell you what is actually present.

8.3.7.1. scipy.signal

sp.signal.sosfilt(sos, x) filters a 1-D signal using cascaded second-order sections (a numerically robust way to apply IIR filters). sos is a sequence of length-6 sections:

from ulab import numpy as np
from ulab import scipy as sp

x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
sos = [[1, 2, 3, 1, 5, 6],
       [1, 2, 3, 1, 5, 6]]
y = sp.signal.sosfilt(sos, x)

You can pass an initial state zi= (shape (n_sections, 2)) to continue filtering across separate buffers. When zi is given, the function returns (y, zf) – the filtered output and the final state – so you can pipe zf into the next call’s zi:

y0, zf0 = sp.signal.sosfilt(sos, buffer0, zi=zi)
y1, zf1 = sp.signal.sosfilt(sos, buffer1, zi=zf0)
# ...

This is the standard way to implement a streaming filter on buffered data.

8.3.7.2. scipy.linalg

Two solver functions:

  • sp.linalg.solve_triangular(a, b, lower=False) – solve a @ x = b with the assumption that a is triangular. The lower keyword picks which triangle to use. Cheap and numerically well-behaved when a is genuinely triangular:

    A = np.array([[3, 0, 0, 0],
              [2, 1, 0, 0],
              [1, 0, 1, 0],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
x = sp.linalg.solve_triangular(A, b, lower=True)

  • sp.linalg.cho_solve(L, b) – given a Cholesky factor L (the kind returned by np.linalg.cholesky), solve A @ x = b where A = L @ L.T. Faster and more accurate than np.linalg.inv when you only want x:

    L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)

8.3.7.3. scipy.optimize

ulab exposes three classic 1-D solvers. Because they have to call back into your Python function on every iteration, the speedup over a pure-Python implementation is modest – around 2x in practice.

8.3.7.3.1. bisect

sp.optimize.bisect(f, a, b, xtol=2.4e-7, maxiter=100) finds a root of a 1-D function by repeatedly halving an interval where the function changes sign:

def f(x):
    return x*x - 1

sp.optimize.bisect(f, 0, 4)
# ~1.0
sp.optimize.bisect(f, 0, 4, maxiter=8)
sp.optimize.bisect(f, 0, 4, xtol=0.1)

8.3.7.3.2. newton

sp.optimize.newton(f, x0, tol=..., rtol=..., maxiter=...) finds a root using a secant / Newton-Raphson iteration starting from x0:

def f(x):
    return x*x*x - 2.0

sp.optimize.newton(f, 3., tol=0.001, rtol=0.01)
# ~1.260...

8.3.7.3.3. fmin

sp.optimize.fmin(f, x0, xatol=..., fatol=..., maxiter=...) finds a local minimum using the downhill-simplex method:

def f(x):
    return (x-1)**2 - 1

sp.optimize.fmin(f, 3.0)
# ~1.0

8.3.7.4. scipy.integrate

Numerical integration, also called quadrature. Four algorithms are available:

  • sp.integrate.quad(f, a, b, order=5, eps=...) – Adaptive Gauss-Kronrod (G10, K21). General-purpose; this is the right default for smooth, well-behaved integrands. Returns a tuple of (value, error_estimate).

  • sp.integrate.romberg(f, a, b, steps=100, eps=...) – classical Romberg / Newton-Cotes method. Returns a single float. Deprecated upstream in favour of quad; included for compatibility.

  • sp.integrate.simpson(f, a, b, steps=100, eps=...) – adaptive Simpson’s rule. Returns a single float.

  • sp.integrate.tanhsinh(f, a, b, levels=6, eps=...) – tanh-sinh (double-exponential) quadrature. The right tool when the integrand has endpoint singularities or one of the limits is infinite. Returns a tuple of (value, error_estimate).

Example – the Gaussian integral:

from math import exp, pi, sqrt
from ulab import numpy as np
from ulab import scipy as sp

f = lambda x: exp(-x*x)
value, err = sp.integrate.tanhsinh(f, -np.inf, np.inf)
print('approx:', value, '   exact:', sqrt(pi))

Note: numerical integration is most accurate when MicroPython is built with double-precision floats (float64). With float32 the routines still work but the achievable tolerance is lower.

8.3.7.5. scipy.special

The special module provides a few statistical / probability functions that work like universal functions: they accept a scalar, an iterable, or an ndarray and return a float ndarray:

from ulab import numpy as np
from ulab import scipy as sp

x = np.linspace(0, 4, num=8)

sp.special.erf(x)         # error function
sp.special.erfc(x)        # complementary error function
sp.special.gamma(x + 1)   # gamma function
sp.special.gammaln(x + 1) # log-gamma function

8.3.7.6. API reference