pythonforbio.
[WASM idle]
NumPy05/18

Fancy Indexing

Starting sandbox…
Beginnerlesson2 tests

Fancy Indexing

In one line: index with an array of integers to gather elements in any order, with repeats — always returns a copy.

The basics

a = np.array([10, 20, 30, 40, 50])

a[[0, 3, 1]]        # array([10, 40, 20])  — arbitrary order
a[[2, 2, 2]]        # array([30, 30, 30])  — repeats allowed
a[np.array([0, -1])] # first and last

Contrast with slicing, which can only express regular strides. Fancy indexing can express any gather, which is why it must copy.

The shape rule

The result has the shape of the index array, not the source array.

a = np.arange(10)
idx = np.array([[0, 1], [2, 3]])    # shape (2, 2)
a[idx].shape                         # (2, 2)

This is what makes fancy indexing a lookup-table operation:

# translate a numeric sequence encoding back to letters
BASES = np.array(["A", "C", "G", "T"])
codes = np.array([0, 3, 3, 2, 1])
BASES[codes]        # array(['A', 'T', 'T', 'G', 'C'])

Multi-dimensional: pairing, not cross-product

This is the trap.

m = np.arange(12).reshape(3, 4)

m[[0, 2], [1, 3]]              # array([1, 11])   ← PAIRS (0,1) and (2,3)
m[np.ix_([0, 2], [1, 3])]      # 2×2 submatrix    ← cross-product

Two index arrays are broadcast against each other and paired elementwise. If you wanted the submatrix formed by rows {0,2} and columns {1,3}, use np.ix_, or index in two steps (m0,2[:, [1,3]], which copies twice).

You can also do it manually with broadcasting:

m[np.array([0, 2])[:, None], np.array([1, 3])[None, :]]   # (2, 1) × (1, 2) → (2, 2)

which is exactly what np.ix_ builds for you. See Broadcasting.

Sorting and ordering

The most common real use of fancy indexing is applying a sort order:

order = np.argsort(pvalues)          # indices that would sort
sorted_p = pvalues[order]
sorted_genes = gene_names[order]     # ← same permutation, labels stay aligned

top10 = np.argsort(pvalues)[:10]
np.argpartition(pvalues, 10)[:10]    # O(n), unsorted top-10 — much faster for large n

np.argmax(x); np.argmin(x)           # single best

argpartition is the one people do not know about. For "give me the 100 most significant genes out of 20,000", it avoids the full sort.

Assignment through fancy indices

a[[0, 2, 4]] = 0
a[[0, 0, 0]] += 1     # ← only adds 1 once, not three times!
np.add.at(a, [0, 0, 0], 1)   # unbuffered — adds three times

The duplicate-index case matters for accumulation. a[idx] += 1 reads, adds, and writes as separate buffered steps, so repeated indices overwrite each other. np.add.at (or np.bincount) does the accumulation properly.

# per-base coverage from a list of read start positions
coverage = np.zeros(chrom_len, dtype=np.int32)
np.add.at(coverage, read_starts, 1)          # correct
# or, faster:
coverage = np.bincount(read_starts, minlength=chrom_len)

take, put and choose

np.take(a, idx)                  # ~ a[idx], but with axis= and mode= options
np.take(m, [0, 2], axis=0)       # explicit axis
np.put(a, idx, values)           # flat scatter
np.take_along_axis(a, idx, axis=1)   # pair argsort with the values it sorted

take_along_axis is the right partner for argsort(axis=...):

order = np.argsort(expr, axis=1)                     # per gene, rank the samples
ranked = np.take_along_axis(expr, order, axis=1)

Bioinformatics examples

# reorder samples to match a metadata table
sample_order = np.array([sample_ids.index(s) for s in metadata["sample"]])
expr_ordered = expr[:, sample_order]

# bootstrap resampling of samples (with replacement)
rng = np.random.default_rng(0)
boot_idx = rng.integers(0, n_samples, size=n_samples)
expr_boot = expr[:, boot_idx]

# one-hot decode
"".join(BASES[onehot.argmax(axis=1)])

# codon table lookup by numeric codon index
aa = CODON_TABLE[codon_indices]

# counting variants per chromosome bin
np.bincount(positions // bin_size, minlength=n_bins)

Fancy vs boolean

Boolean mask Integer array
Length must match the axis any length
Order preserves original order arbitrary
Repeats impossible allowed
Use for filtering reordering, lookup, resampling

np.flatnonzero(mask) converts a mask to indices when you need to combine both.

Common mistakes

  • Expecting m[rows, cols] to be a submatrix. It pairs. Use np.ix_.
  • a[dup_idx] += 1 for accumulation. Use np.add.at or np.bincount.
  • Forgetting fancy indexing copies. On a large matrix, expr[big_idx_array] allocates a whole new matrix — watch memory.
  • Applying a sort order to the data but not the labels. Same failure mode as in Boolean Masking.
  • Using a float array as an index. IndexError. Cast with .astype(int).
  • Full argsort when you only need the top k. Use argpartition.

See also

Indexing and Slicing · Boolean Masking · Views vs Copies · Broadcasting · Sorting and Ranking · NumPy Random Generator

Test cases · 2

#viainputexpected stdout
1stdinTP53 BRCA1 EGFR KRAS 120 45 300 300 2['EGFR', 'KRAS']
2stdinA B C 1 2 3 3['C', 'B', 'A']

Hints · 2

01Hint

An array of integers used as an index selects those positions in order.

02Hint

np.argsort returns the indices that would sort an array.

starter

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

BroadcastingIndexing and Slicing