The array and array_ufunc Protocols

The __array__ method tells NumPy how to convert your object to an ndarray. NumPy calls this whenever it needs to operate on your object:

class PhysicalQuantity:
    """A value with units that converts to raw NumPy array."""

    def __init__(self, value: np.ndarray, unit: str = "V"):
        self._value = np.asarray(value)
        self.unit = unit

    def __array__(self, dtype=None, copy=None) -> np.ndarray:
        if dtype is not None:
            return self._value.astype(dtype)
        return self._value

    def __repr__(self) -> str:
        return f"PhysicalQuantity({self._value}, unit='{self.unit}')"

Now np.array(PhysicalQuantity([1, 2, 3])) works, and so does np.mean(PhysicalQuantity([1, 2, 3])).

__array_ufunc__ intercepts NumPy universal function (ufunc) calls on custom objects. It gives your class full control over how operations like np.add, np.multiply, and np.sin are handled:

class PhysicalQuantity:
    # ... (previous code) ...

    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
        """Intercept NumPy ufunc calls to preserve units."""
        # Convert all PhysicalQuantity inputs to raw arrays
        raw_inputs = []
        units = set()
        for inp in inputs:
            if isinstance(inp, PhysicalQuantity):
                raw_inputs.append(inp._value)
                units.add(inp.unit)
            else:
                raw_inputs.append(inp)

        # Apply the ufunc to raw arrays
        result = getattr(ufunc, method)(*raw_inputs, **kwargs)

        # Wrap result back in PhysicalQuantity
        if len(units) == 1:
            return PhysicalQuantity(result, unit=units.pop())
        return result  # mixed units: return raw

The signature is always (self, ufunc, method, *inputs, **kwargs) where method is typically '__call__', 'reduce', 'accumulate', or 'outer'.

While __array_ufunc__ handles element-wise operations (ufuncs), __array_function__ (NEP 18) intercepts all NumPy functions including np.concatenate, np.linalg.norm, np.fft.fft, etc.:

import functools

HANDLED_FUNCTIONS = {}

def implements(np_function):
    """Register an implementation for a NumPy function."""
    def decorator(func):
        HANDLED_FUNCTIONS[np_function] = func
        return func
    return decorator

class PhysicalQuantity:
    def __array_function__(self, func, types, args, kwargs):
        if func not in HANDLED_FUNCTIONS:
            return NotImplemented
        return HANDLED_FUNCTIONS[func](*args, **kwargs)

@implements(np.concatenate)
def concatenate(arrays, axis=0):
    raw = [np.asarray(a) for a in arrays]
    units = {a.unit for a in arrays if isinstance(a, PhysicalQuantity)}
    result = np.concatenate(raw, axis=axis)
    if len(units) == 1:
        return PhysicalQuantity(result, unit=units.pop())
    return result

This is how libraries like pint and astropy.units provide transparent NumPy integration.