pythonforbio.
[WASM idle]
NumPy18/18

ufuncs

Starting sandbox…
Beginnerlesson2 tests

ufuncs

In one line: universal functions — compiled, elementwise, broadcasting-aware operations that are the reason NumPy is fast.

Every +, *, np.log2, np.exp and np.maximum you write is a ufunc. Understanding their shared machinery gives you several free capabilities most people never use.

What counts as a ufunc

np.add, np.subtract, np.multiply, np.divide, np.power, np.mod
np.exp, np.log, np.log2, np.log10, np.log1p, np.sqrt, np.abs
np.sin, np.cos, np.tan
np.maximum, np.minimum          # elementwise pairwise — NOT max/min reductions
np.greater, np.equal, np.logical_and
np.isnan, np.isinf, np.isfinite
np.floor, np.ceil, np.round, np.clip

Operators dispatch to them: a + b is np.add(a, b).

np.maximum vs np.max is a classic confusion. np.maximum(a, b) compares two arrays elementwise; np.max(a) reduces one array to its largest value. Same for minimum/min.

The shared keyword arguments

Every ufunc accepts these, and they are underused:

np.add(a, b, out=result)              # write into an existing buffer, no allocation
np.log2(x, where=x > 0)               # only compute where the condition holds
np.divide(1, x, out=np.zeros_like(x), where=x != 0)   # safe reciprocal
np.add(a, b, dtype=np.float64)        # force accumulation precision

The out= + where= combination is the correct way to avoid divide-by-zero and log-of-zero warnings. The common alternative, np.where(x > 0, np.log2(x), 0), still computes log2(0) and emits the warning, because both branches are evaluated eagerly.

out= also matters in memory-tight loops:

np.log2(big, out=big)      # transform in place, no second 8 GB allocation
big += 1                    # likewise

The four ufunc methods

Binary ufuncs carry four extra methods that give you reductions for free:

np.add.reduce(a)                    # ≡ a.sum()
np.add.accumulate(a)                # ≡ a.cumsum()
np.add.reduceat(a, indices)         # segmented reduction — grouped sums!
np.add.outer(a, b)                  # every pair — same as a[:,None] + b[None,:]
np.add.at(a, idx, values)           # unbuffered in-place scatter — handles duplicates

np.multiply.reduce(a)               # ≡ a.prod()
np.maximum.accumulate(a)            # running maximum
np.logical_or.reduce(masks)         # OR together a stack of masks

np.add.reduceat is a genuine grouped-sum primitive:

# sum coverage within exons, given exon start indices in a sorted array
starts = np.array([0, 150, 400])
np.add.reduceat(coverage, starts)     # one total per segment

np.maximum.accumulate is how you compute a running maximum in one pass, useful for step-down multiple-testing procedures.

Type promotion and casting

np.add(np.int8(100), np.int8(100))     # 200 → overflows to -56 in int8
np.int32([1]) + np.float32([1])        # → float64 (int32 needs more precision than float32 has)
np.float32([1]) + 1.0                  # → float32 under NEP 50 (NumPy 2)

The NumPy 2 promotion change (NEP 50) means Python scalars no longer upcast arrays. See NumPy dtypes.

Speed, concretely

x = np.random.default_rng(0).random(1_000_000)

[math.log2(v) for v in x]     # ~250 ms  — 1M interpreter round-trips
np.log2(x)                     # ~8 ms    — one C loop, SIMD-vectorised

Roughly 30× here; for simpler operations the gap widens to 100×+. The reason is not that C is faster than Python per operation — it is that the loop, the type dispatch and the boxing all happen once instead of a million times. See Vectorization.

Writing your own

If you have a scalar function you cannot express in NumPy operations:

@np.vectorize            # ← convenience, NOT speed. It is a Python loop.
def f(x): ...

from numba import vectorize
@vectorize(["float64(float64)"], nopython=True)
def f(x): ...            # ← this one is genuinely compiled

np.vectorize is not a performance tool. Its docstring says so. It exists to give a scalar function ufunc-like broadcasting semantics. If you need speed, use numba, or restructure the problem into existing ufuncs.

Bioinformatics examples

# Phred quality ↔ error probability
p_err = 10 ** (-quals / 10)
quals = -10 * np.log10(p_err)

# log transform with a pseudocount (the standard RNA-seq move)
np.log2(counts + 1)

# safe log with an explicit floor
np.log2(counts, out=np.full_like(counts, -np.inf), where=counts > 0)

# clip outliers before plotting
np.clip(log2fc, -5, 5)

# elementwise max against a floor
np.maximum(pvalues, 1e-300)      # avoid log10(0) = -inf in a volcano plot

# GC content across a one-hot matrix
onehot[:, [1, 2]].sum() / onehot.sum()

Common mistakes

  • np.max when you meant np.maximum (or vice versa). One reduces, one pairs.
  • Expecting np.where to short-circuit. Both branches are always computed. Use where= on the ufunc.
  • np.vectorize for speed. It is a loop.
  • Ignoring RuntimeWarnings for invalid value / divide by zero. They mean NaN or inf entered your data. Use np.errstate to handle deliberately, not warnings.filterwarnings to hide:
    with np.errstate(divide="ignore", invalid="ignore"):
        ratio = a / b
    
  • out= with an incompatible dtype. It silently casts, which can truncate.
  • Forgetting integer overflow in np.add.reduce on int32 counts. Pass dtype=np.int64.

See also

Vectorization · Broadcasting · Axes and Reductions · NumPy dtypes · Missing Data in NumPy

Test cases · 2

#viainputexpected stdout
1stdin0 1 3 7 15[0.0, 1.0, 2.0, 3.0, 4.0]
2stdin31 63[5.0, 6.0]

Hints · 2

01Hint

np.log2 is a ufunc: it applies element-wise with no Python loop.

02Hint

Add the pseudocount before taking the log, not after.

starter

No output yet — run the code to populate this drawer.

ndarrayend of path