ufuncs

Core Principle: Universal functions (ufuncs) are compiled, vectorized C routines that execute elementwise computations with built-in broadcasting and type resolution.

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

(<ufunc 'floor'>, <ufunc 'ceil'>, <function round at 0x10623af30>, <function clip at 0x1062397f0>)

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 across contiguous interval boundaries
np.add.outer(a, b)                  # outer product across pairwise combinations
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

array([2.], dtype=float32)

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

Speed, concretely

array([-0.6507215 , -1.89010879, -4.60916421, ..., -2.91210759,
       -2.06869819, -1.03919183], shape=(1000000,))

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

Exercise

Write log_ratios(a, b), which computes log2(a / b) elementwise but only where both operands are strictly positive. Use the out and where keywords of np.divide and np.log2 so the undefined positions are never evaluated at all; an implementation that divides first and masks afterwards raises warnings and produces infinities.

Return five lines. The first is ratios=<values joined by commas>, each defined value rounded to three decimals and each undefined position written as NA. The second is defined=<n>/<total>. The third is sum=<value> and the fourth product=<value>, computed over the defined values with np.add.reduce and np.multiply.reduce and rounded to three decimals; when nothing is defined these must be the identities of the two operations, which are not the same number. The fifth is running=<list>, the cumulative totals from np.add.accumulate over the defined values, rounded to three decimals, and the empty list when there are none.