9.16. Numerical extras

The pages so far covered the everyday surface – construction, indexing, math, broadcasting, reductions, linear algebra, FFT, filtering. This page catalogues the remaining numerical tools that did not fit one of those headings: interpolation, polynomial fitting, integration, convolution, root finding, optimisation, special functions, random numbers.

Almost every function on this page accepts – and most require – ndarray inputs and returns either a float scalar or a float ndarray. The imports:

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

9.16.1. Interpolation

interp(x, xp, fp)() does one-dimensional linear interpolation. xp is a monotonically increasing 1-D array of independent values; fp is the matching dependent values; x is where the interpolant should be evaluated:

xp = np.array([0.0, 1.0, 2.0, 3.0])
fp = np.array([10.0, 20.0, 25.0, 30.0])
x  = np.array([0.5, 1.5, 2.5])

np.interp(x, xp, fp)
# array([15.0, 22.5, 27.5])

Outside the [xp[0], xp[-1]] range the result is clipped to fp[0] and fp[-1] respectively; the left= and right= keywords override those endpoints.

Used on the camera to remap a calibration table to arbitrary sample positions – a temperature-to-voltage table from a thermistor, a non-linear pixel response curve, a per-channel gamma lookup. One library call, one pass over the inputs, no Python loop.

9.16.2. Polynomial fitting and evaluation

polyfit(x, y, deg)() fits a polynomial of degree deg to the data points (x, y) by least squares and returns the coefficients (highest degree first):

x = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
y = np.array([0.0, 1.1, 3.9, 9.1, 15.8])

coeffs = np.polyfit(x, y, 2)
# array([1.0, 0.0, 0.0])  (approximately)

When x is omitted, range(len(y)) is used – useful for fitting a quick polynomial to a regularly sampled buffer with no associated x-axis:

np.polyfit(y, 2)

polyval(p, x)() evaluates the polynomial whose coefficients are p at x. The input x can be a scalar (returns a float) or an ndarray (returns an ndarray):

p = np.polyfit(x, y, 2)
fitted = np.polyval(p, x)

9.16.3. Convolution

convolve(a, v)() returns the full-length discrete linear convolution of two 1-D arrays. Only 'full' mode is implemented; the output length is len(a) + len(v) - 1. Slice the result for the same effect as the 'same' and 'valid' modes that the desktop numpy offers:

a = np.array([1.0, 2.0, 3.0])
v = np.array([0.5, 0.5])

np.convolve(a, v)
# array([0.5, 1.5, 2.5, 1.5])

Useful for short FIR filters and smoothing kernels (box, triangle, gaussian) where setting up an SOS chain is overkill.

9.16.4. Trapezoidal integration

0)() integrates a sampled function by the composite trapezoidal rule:

x = np.linspace(0, np.pi, num=128)
y = np.sin(x)
np.trapz(y, x)             # ~ 2.0

Pass dx= when the sample spacing is uniform and only the step matters; pass x= when the samples are not evenly spaced. The right call for integrating already- captured sensor data, where the analytic form is not available.

9.16.5. Numerical integration of a callable

When the integrand is a Python function rather than a buffer of samples, scipy.integrate exposes four quadrature algorithms:

  • )() – adaptive Gauss-Kronrod. The right default for smooth integrands. Returns (value, error).

  • )() – classical Romberg / Newton-Cotes. Returns a single float. Deprecated upstream; included for compatibility.

  • )() – adaptive Simpson’s rule. Returns a single float.

  • )() – double-exponential quadrature. Use when the integrand has endpoint singularities or an infinite limit. Returns (value, error).

The Gaussian integral as a worked example:

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))

Output:

approx: 1.7724538...   exact: 1.7724538...

Numerical integration is most accurate on cams whose firmware uses double-precision floats. With single precision the routines still work, but the achievable tolerance is lower.

9.16.6. Root finding and minimisation

scipy.optimize covers three classic single-variable solvers. Each iteration calls back into the user-supplied Python function, so the speedup over a pure-Python solver is modest (roughly 2x); the convenience is in not having to write the solver.

  • )() – find a root of f on [a, b] by halving the interval. f(a) and f(b) must have opposite signs:

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

sp.optimize.bisect(f, 0, 4)        # ~1.0

  • )() – find a root using secant / Newton-Raphson iteration:

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

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

  • )() – find a local minimum using the downhill-simplex (Nelder-Mead) method:

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

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

9.16.7. Special functions

scipy.special exposes a handful of statistical and probability functions that behave like universal functions – they accept a scalar, an iterable, or an ndarray and return a float ndarray:

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

The error function and its complement appear in the CDF of a Gaussian; the gamma functions show up in beta / chi-squared / student-t calculations.

9.16.8. Random numbers

numpy.random provides a Generator class that draws samples from common distributions. The generator is stateful: each call advances its internal state, so consecutive calls return independent samples:

from ulab import numpy as np

rng = np.random.Generator(seed=42)

rng.random(size=5)             # 5 uniform [0.0, 1.0) samples
rng.uniform(low=-1.0, high=1.0, size=10)
rng.normal(loc=0.0, scale=1.0, size=(2, 4))

The output dtype is always float. size= accepts an integer (1-D output) or a tuple (n-D output); when omitted, a single Python float is returned.

The generator is suitable for simulation, dithering, synthetic test data, and any other application where cryptographic strength is not required. It is not suitable for keys or tokens.

9.16.9. What is not covered here

A handful of routines exposed by the reference – the .npy I/O helpers (load(), save(), loadtxt(), savetxt()), set_printoptions(), the experimental curve_fit() stub – are straightforward extensions of patterns already shown on this page. See numpy — numpy-compatible array operations and scipy — subset of scipy via ulab for the complete argument-level reference.

9.16.10. Firmware availability

Whether each submodule is actually present depends on how the firmware on the cam was built. Calling a function the firmware does not include raises AttributeError. dir(sp), dir(sp.optimize), dir(np.random) and friends report what is available on the cam being targeted.