# NumPy

Source: https://codewiki.com/datascience/numpy/

> - **what**: NumPy represents homogeneous multidimensional data with `ndarray` and applies numerical operations across whole arrays.
> - **when**: Use it when data forms a regular numeric grid and you need slicing, aggregation, linear algebra, or interchange with Python's scientific libraries.
> - **how**: Define what each axis and `dtype` mean, then compose array operations; check the result shape and memory ownership around broadcasting, slicing, and conversion.

## What it is and why it exists

NumPy is Python's numerical array library. Its core object is the N-dimensional array (`ndarray`): homogeneous elements plus metadata that describes dimensions, data type, and memory layout. A two-dimensional array might mean "samples × features," while a three-dimensional array might mean "image height × width × channels," but NumPy doesn't record those domain meanings for you.

Python lists are good at holding heterogeneous objects and changing collections, but they don't have one set of array arithmetic rules. `[1, 2] * 2` repeats a list, whereas multiplying a NumPy array by `2` computes element by element. NumPy gives regular numeric data a shared array protocol, so Pandas, SciPy, Matplotlib, and many machine learning tools can exchange data without a Python loop for every operation.

NumPy fits regular, dense, homogeneous numeric grids. Tables with column names and mixed types usually fit Pandas or Polars better; highly irregular object collections may still belong in Python containers. That boundary matters: don't force irregular data into a `dtype=object` array merely to call the result "vectorized."

The examples here were run with Python 3.14.3 and NumPy 2.5.2. Their displayed values use explicit rounding or list conversion, so they don't depend on array truncation rules or terminal width.

## How it works

### Three contracts in one array

Every array has an array shape, a data type (`dtype`), and a data buffer. `shape` is a tuple containing the length of each axis; `ndim` is the number of axes; `size` is the total element count. The `dtype` tells NumPy how to interpret each element, such as `int64`, `float32`, or `bool`.

An array variable also has an implicit domain contract: what each axis means. Two arrays with shape `(32, 10)` could represent 32 samples with 10 features or 32 time points from 10 sensors. NumPy checks whether lengths can participate in an operation, not whether the business meanings of the axes have been reversed.

An array stride is the number of bytes crossed when moving one position along an axis. Slicing and transposition can change only `shape` and `strides` while continuing to reference the same buffer. Two arrays that look different can therefore share data, which is why modifying a slice sometimes changes its source.

### Creation fixes the initial type

`np.array()` builds from existing Python data; `np.zeros()`, `np.ones()`, and `np.full()` create initialized regular arrays; `np.arange()` and `np.linspace()` create numeric sequences. Specify `dtype` at business boundaries, especially when integer width and floating precision matter or when input may be empty. This prevents one batch of sample values from deciding the type by accident.

`np.empty()` allocates space without initializing its values. It suits internal code that will overwrite every element, not an "empty values" array. If any position could be read before assignment, use a creation function with a defined initial value.

`reshape()` requires the element count to stay unchanged. `-1` asks NumPy to infer one axis length, but a shape-valid statement can still express the wrong data layout. Name the meaning of each axis before and after a reshape instead of checking only whether the element count divides evenly.

Common shape operations differ in useful ways:

| Goal | Expression | Result |
| --- | --- | --- |
| Change dimension lengths | `a.reshape(rows, cols)` | Keeps the element count and may share memory |
| Swap all axes | `a.T` | Swaps rows and columns in 2D and is usually a view |
| Name an axis order | `a.transpose(order)` | `order` maps output axes to input axes |
| Move one axis | `np.moveaxis(a, source, destination)` | Puts a channel or similar axis in an explicit position |
| Insert a length-one axis | `np.expand_dims(a, axis)` | Often prepares an array for broadcasting |
| Remove a length-one axis | `np.squeeze(a, axis=axis)` | Naming `axis` avoids deleting other axes |
| Combine along a new axis | `np.stack(arrays, axis=axis)` | Requires every input to have the same shape |
| Join along an existing axis | `np.concatenate(arrays, axis=axis)` | Requires other axis lengths to match |

### Moving and joining axes

Transposition doesn't mean "swap rows and columns for every array." `.T` swaps the two axes of a 2D array, but it reverses the complete axis order of a 3D array. For batch, channel, and spatial dimensions, `transpose()` or `moveaxis()` states the intended order more clearly.

`stack()` and `concatenate()` solve different problems. `stack()` creates a new axis, so ten images shaped `(64, 64)` become `(10, 64, 64)`; `concatenate()` extends an existing axis. Generated code often adds or loses one batch axis at this boundary.

Calling `squeeze()` without an `axis` removes every length-one axis. That seems convenient in sample data, but it can make a return value lose a dimension whenever the batch size happens to be one. Public APIs should name the axis to remove or keep a stable batch dimension.

`np.newaxis` is an alias for `None` that inserts a length-one axis in an index. `values[:, np.newaxis]` turns `(n,)` into `(n, 1)`, while `values[np.newaxis, :]` produces `(1, n)`. The elements are identical, but the broadcasting direction is completely different.

Before joining arrays, validate lengths outside the target axis. An error can report that dimensions conflict, but it can't know that one input is channel-first while another is channel-last. Axis names in the interface matter more than repairing one `ValueError`.

### Vectorization, universal functions, and reductions

Vectorization expresses a computation as operations on complete arrays. For example, `prices * quantities` means element-wise multiplication, while `np.sqrt(values)` applies one universal function to every element. Vectorization first changes the expression of the problem: NumPy handles the loop, and shape and type rules live in the array operation.

Reductions such as `sum()`, `mean()`, and `max()` use `axis` to name the axis they remove. For a "samples × features" matrix, `axis=0` removes the sample axis and leaves one result per feature; `axis=1` removes the feature axis and leaves one result per sample. `keepdims=True` retains a length-one axis and often makes later broadcasting easier to read.

Ordinary arithmetic is element-wise. `*` between two matrices is not matrix multiplication; use `@` or `np.matmul()` for linear algebra multiplication. Use `np.linalg.solve()` for a square linear system rather than explicitly forming an inverse and multiplying by it.

### Broadcasting starts at the trailing axes

Broadcasting lets arrays with different shapes participate in element-wise operations. NumPy compares axes from the right: a pair is compatible when the lengths are equal or one length is `1`; a missing leading axis behaves as length `1`. Otherwise, the operation raises a shape-compatibility error.

For example, a sample matrix shaped `(4, 3)` can subtract feature means shaped `(3,)` because both trailing axes have length `3`. Per-sample means shaped `(4,)` can't be subtracted directly because trailing lengths `3` and `4` conflict. Retaining those means as `(4, 1)` gives them a length-one trailing axis that broadcasts across features.

Broadcasting describes a conceptual expansion; it doesn't mean NumPy first creates repeated input data. The arithmetic result will usually still allocate its own array. Predict the result shape before writing the expression. "It didn't raise" proves only that the shapes are compatible, not that the axes mean the right things.

### Indexing decides sharing or copying

Basic slicing, such as `array[1:4]` or `matrix[:, ::2]`, usually returns an array view that shares data. Element changes through that view are visible in the source. `view.base` helps exploration, but the base of a chain of views isn't guaranteed to be the array in your hand, so use `np.shares_memory()` to test whether two arrays overlap.

Integer-array and Boolean indexing are advanced indexing and return copies. You can modify `array[[0, 2]]` and `array[array > 0]` independently afterward, but the copy allocates a new buffer. When basic and advanced indexing are mixed, the result shape can also depart from a simple axis-by-axis intuition, so confirm it with a small input.

Assignment has a separate rule: `array[mask] = 0` writes into the source because the index is the assignment target. In contrast, assigning `selected = array[mask]` and then changing `selected` changes only that copy. A review must follow the complete data flow; one index expression isn't enough to tell where the final write lands.

### Masks, conditions, and order

Comparing arrays produces Boolean arrays. Combine conditions with the element-wise operators `&`, `|`, and `~`, usually with parentheses around each comparison. Python's `and` and `or` require one truth value for the complete array and therefore raise an ambiguity error.

Boolean indexing retains only matching elements and often flattens the original structure. `np.where(condition, left, right)` instead chooses element by element between two equally shaped or broadcastable values and keeps the broadcasted shape. The right choice depends on whether downstream code still needs the original axes.

Missing data has more than one representation. Floating arrays can store `NaN` directly, while integer arrays cannot; creating an array from integers and `NaN` will generally promote it to a floating type. If missingness is mixed with zero, an empty string, or a sentinel, normalize that meaning before a reduction.

`np.argsort()` returns indexes that can reorder an array, not the sorted values themselves. To reorder labels and features together, compute one order and apply it to every associated array. Sorting only one column breaks the correspondence within each sample.

Mask and order-index lengths are part of the array contract too. A stale mask produced before one filter can't safely be reused after that filter. A pipeline should carry a mask alongside its matching data or recompute it at each stage.

### Errors and numerical warnings

Incompatible shapes usually raise `ValueError`. Read the input shapes in the exception and compare them from the trailing axes before inserting a blind `reshape()`. A shape patch that makes code run doesn't necessarily preserve the data's meaning.

Unsafe in-place casts fail. For example, a floating division result can't be written back into an integer array under the `same_kind` casting rule. The usual fix is to create a result with the required floating type, not force a conversion that truncates information.

Floating divide-by-zero, invalid operations, and overflow may emit `RuntimeWarning` by default while producing `inf` or `NaN`. Those values can continue into later calculations if the caller ignores warnings, so "the function didn't raise" isn't a sufficient check.

`np.errstate()` changes floating-error handling inside one explicit block. Tests can temporarily promote relevant warnings to exceptions on critical numerical paths. An algorithm that expects special values can handle them locally and validate the result immediately.

`np.seterr()` changes NumPy's error settings for the current thread and shouldn't be changed casually inside an ordinary library function. A local context is easier to reason about than hidden environmental state and doesn't make the caller restore settings afterward.

`np.isfinite()` finds `NaN` and both infinities together. Input validation and important intermediate results can use it, but the domain contract still decides whether those values are rejected, replaced, or preserved.

Warning policy doesn't replace boundary tests. Cover division by zero, very large and small values, and all-missing slices separately because their promotion, warnings, and result shapes can differ.

## Examples

### Read shape, type, and axes

Rows in this matrix are days and columns are sensors. `axis=1` removes the sensor axis to produce daily means; `axis=0` removes the day axis to produce a peak for each sensor.

<!-- quick -->

```python
# file: array_axes.py
import numpy as np

temperatures = np.array(
    [[18.5, 20.0, 19.5], [21.0, 23.5, 22.0]],
    dtype=np.float64,
)

print(f"shape={temperatures.shape}")
print(f"ndim={temperatures.ndim}, dtype={temperatures.dtype}")
print("daily mean:", np.round(temperatures.mean(axis=1), 2).tolist())
print("sensor peak:", temperatures.max(axis=0).tolist())
```

```text
shape=(2, 3)
ndim=2, dtype=float64
daily mean: [19.33, 22.17]
sensor peak: [21.0, 23.5, 22.0]
```

<!-- /quick -->

The shape `(2, 3)` doesn't say which axis is the day. Variable names, nearby comments, and interface documentation must preserve that meaning. After a reduction, check that the result shape matches what its consumer expects.

### Standardize features with broadcasting

This example treats each column as one feature. `keepdims=True` keeps the mean and standard deviation shaped `(1, 3)`, so subtraction and division clearly broadcast across samples. The third column is constant, and the code replaces its zero standard deviation with `1.0` to avoid division by zero.

```python
# file: standardize_features.py
import numpy as np

samples = np.array(
    [[1.0, 10.0, 7.0], [2.0, 20.0, 7.0], [3.0, 30.0, 7.0]]
)

means = samples.mean(axis=0, keepdims=True)
scales = samples.std(axis=0, keepdims=True)
safe_scales = np.where(scales == 0, 1.0, scales)
standardized = (samples - means) / safe_scales

np.set_printoptions(precision=2, suppress=True)
print("mean shape:", means.shape)
print(standardized)
print("constant feature:", standardized[:, 2].tolist())
```

```text
mean shape: (1, 3)
[[-1.22 -1.22  0.  ]
 [ 0.    0.    0.  ]
 [ 1.22  1.22  0.  ]]
constant feature: [0.0, 0.0, 0.0]
```

Replacing a zero scale with `1.0` is an explicit policy in this example, not a universal rule. Another system might drop constant features or reject the input. Numerical code needs to put that domain choice in its interface and tests.

### Distinguish a view from a copy

The basic slice `scores[:, 1:]` shares data with the source, while Boolean indexing creates a copy. The assignments look similar, but only the first one writes back into `scores`.

```python
# file: view_copy.py
import numpy as np

scores = np.array(
    [[70, 80, 90], [60, 75, 85], [88, 92, 95]],
    dtype=np.int64,
)

last_two = scores[:, 1:]
selected = scores[scores[:, 0] >= 70]

print("slice shares memory:", np.shares_memory(scores, last_two))
print("mask shares memory:", np.shares_memory(scores, selected))

last_two[0, 0] = 999
selected[0, 0] = -1
print("source first row:", scores[0].tolist())
print("selected first row:", selected[0].tolist())
```

```text
slice shares memory: True
mask shares memory: False
source first row: [70, 999, 90]
selected first row: [-1, 80, 90]
```

Call `.copy()` when a modification must be isolated in a new buffer. When a function updates in place, say so in its name or documentation. Making callers guess about ownership lets an apparently local cleaning step corrupt upstream data.

### Isolate random state and solve least squares

`default_rng()` returns an independent random generator rather than resetting module-level random state. This example uses it to produce training and test indexes, then fits a one-feature linear model with an intercept through `np.linalg.lstsq()`.

```python
# file: random_lstsq.py
import numpy as np

rng = np.random.default_rng(42)
feature = np.array([0, 1, 2, 3, 4, 5], dtype=np.float64)
target = 1.0 + 2.0 * feature

order = rng.permutation(feature.size)
train_index = order[:4]
test_index = order[4:]

design = np.column_stack(
    [np.ones(train_index.size), feature[train_index]]
)
intercept, slope = np.linalg.lstsq(
    design, target[train_index], rcond=None
)[0]
predicted = intercept + slope * feature[test_index]

print("train index:", train_index.tolist())
print(f"model: y = {intercept:.1f} + {slope:.1f}x")
print("test index:", test_index.tolist())
print("prediction:", np.round(predicted, 6).tolist())
```

```text
train index: [3, 2, 5, 4]
model: y = 1.0 + 2.0x
test index: [1, 0]
prediction: [3.0, 1.0]
```

A fixed seed makes this run repeatable, but a long-lived experiment should also save the actual indexes, NumPy version, and selected bit generator. Don't treat "same seed" as complete data provenance across every version and algorithm. Real model splits must also define constraints such as stratification, grouping, and time order before shuffling.

## Pitfalls

> **Pitfall:** Treating "broadcasting succeeded" as "the axes aligned." Adding shapes `(100, 1)` and `(100,)` produces `(100, 100)`, even if the intent was to add one value per sample.

**Fix:** Assert shapes at boundaries and document what every axis means. Use `keepdims=True`, `np.newaxis`, or an explicit `reshape()` to name the expanded axis, then check the result with a tiny hand-calculated input.

> **Pitfall:** Assuming slices are always independent, or assuming every index returns a view. Basic slices usually share memory, advanced indexes usually copy, and mixed indexing is harder to judge by appearance.

**Fix:** Call `.copy()` when isolation is required and `np.shares_memory()` when overlap must be confirmed. Tests for in-place functions should check both the returned value and whether the input changed according to contract.

> **Pitfall:** Ignoring fixed-width integer overflow or writing a floating result back into an integer array. Array operations follow NumPy's type rules and don't automatically gain the arbitrary precision of Python integers.

**Fix:** Choose a wide enough `dtype` at the input boundary and use `np.result_type()` to inspect the result type of a mixed operation. Check types before division, normalization, and cumulative sums, and test maxima, minima, and negative values.

> **Pitfall:** Treating `np.empty()` as an array filled with empty values or zeros. It can expose arbitrary bit patterns left in that memory block, so an unwritten position has no domain meaning.

**Fix:** Use `empty()` only when you can prove every element is assigned before it is read. Otherwise, use `zeros()`, `full()`, or an explicit initial value that communicates the missing-data policy.

> **Pitfall:** Calling global `np.random.seed()` inside a library function. It resets legacy random state shared with other code in the process and makes tests depend on call order.

**Fix:** Accept a `Generator` parameter or create a local `Generator` with an explicit seed inside the function. When an experiment must replay across processes, save the samples and bit-generator information too.

> **Pitfall:** Comparing floating calculations with `==`, or letting one `NaN` silently propagate through a whole reduction. Binary floating-point rounding and missing-data policy are separate problems.

**Fix:** Use `np.isclose()` or `np.allclose()` with tolerances from the domain. Decide whether `NaN` should be rejected, ignored, or propagated before choosing `mean()`, `nanmean()`, or explicit validation; don't ignore missing values merely to obtain a number.

<!-- deep -->

## Memory layout and data ownership

An `ndarray` separates a data buffer from the metadata that interprets it. That metadata includes `shape`, `dtype`, `strides`, and an offset from the start of the buffer. A view can change those fields without copying elements, so a transpose, reverse slice, or stepped slice can still see the same underlying data.

This design means contiguity isn't guaranteed. A C-contiguous array usually has adjacent elements along its last axis; a Fortran-contiguous array usually has them along its first. A transposed array may satisfy one layout or neither. Read an external library's contract and use `np.ascontiguousarray()` to materialize a contiguous copy only when the interface requires one.

`ravel()` returns a view when it can but may copy; `flatten()` always returns a copy. `reshape()` can produce a view when the layout permits it, may copy otherwise, and can fail for some shape and order combinations. When ownership matters, don't replace checks and tests with a claim that a function "usually doesn't copy."

A broadcasted view may map several logical positions to one element. `np.broadcast_to()` returns a read-only view so one assignment can't ambiguously affect several positions. Arithmetic results are usually new arrays, but an `out=` parameter or in-place operation changes the ownership and casting requirements and should be documented at the interface.

You can audit element-buffer size directly: `size * itemsize` equals the array's `nbytes`. That figure excludes the small Python object metadata and doesn't mean a view owns another buffer of the same size. When investigating peak memory, distinguish a view's logical size, the buffer it owns, and temporary results created by an expression.

## Type promotion and numerical boundaries

`dtype` is part of an operation's behavior, not merely a storage choice. When two inputs participate, NumPy derives a common result type from their data types and the kinds of any Python scalars; `np.result_type()` can inspect that decision without performing the complete calculation. Assignment into an existing array must also satisfy casting rules, so a floating result can't safely be written in place to an integer array.

Fixed-width integers overflow at their boundaries. A reduction may need a wider range than each individual input, so counts, pixels, or small integer features can fit in `int8` while their sum cannot. Choose a type for the final operation's range, not only for the raw values.

Floating types require decisions about range, precision, and special values. `NaN` doesn't equal itself, positive and negative infinity can be found with `np.isfinite()`, and a near-zero denominator magnifies rounding error. Tolerances and missing-data policy come from the problem. NumPy supplies checks but doesn't decide what result is acceptable.

Array APIs also distinguish element-wise semantics from linear algebra semantics. `a * b`, `a @ b`, and `np.dot(a, b)` need not agree for higher-dimensional inputs. Public functions should constrain accepted dimensions and test non-square batches instead of inferring high-dimensional behavior from a two-dimensional example.

Randomness has state and algorithm boundaries too. Accepting a `Generator` as a parameter gives ownership of random state to the caller and lets tests inject a fixed generator. Exact experimental reproduction calls for saving inputs, actual samples, bit generator, and library version, not just one integer seed.

## Array contracts at API boundaries

### Input validation

A public function that accepts arrays should decide whether it takes any "array-like" object or only an `ndarray`. `np.asarray()` can normalize lists and similar inputs and avoid a needless copy when the input is already compatible. After normalization, check `ndim`, `shape`, `dtype`, and finite values so errors refer to one representation.

Don't write only `assert array.shape[1] == 3`. Ordinary assertions can be removed in optimized mode, and a 1D input will raise an indexing error before that condition can explain the problem. Validate type, dimensions, axis lengths, and numeric range in order, with an exception message that tells the caller how to correct the input.

Empty arrays need a separate domain decision. Some element-wise transformations can naturally return an empty result, while a reduction such as `max()` has no answer without an identity. State whether the API rejects empty input, returns an empty array, or uses an explicit `initial` value.

Writeability is also part of the contract. Setting `array.flags.writeable` to `False` can prevent accidental in-place changes, but it isn't a security boundary and doesn't freeze other shared views. When isolation is required, make a copy and control who keeps writeable references.

### Output shape and ownership

Function documentation should state an output shape relative to its input, such as "input `(samples, features)`, output `(features,)`." Saying only "returns the mean" doesn't distinguish a global scalar, per-sample means, and per-feature means. Type annotations don't yet replace this concrete axis contract.

Returning a view can avoid a copy, but it exposes the input's lifetime and mutability to the caller. Returning a copy isolates the result but allocates a new buffer. Either choice can be sound; it needs to be stable, visible, and covered by tests that modify both input and output.

If a function accepts an `out=` parameter, validate the target's shape, writeability, and casting safety. Overlapping input and output can also change the computation. Unless the underlying operation explicitly supports that overlap, don't assume an in-place write is equivalent to computing the whole result before assignment.

Scalar return values deserve attention as well. Some NumPy reductions return a NumPy scalar such as `np.float64`; it participates in most Python arithmetic, but a serializer may require a native `float` or `int`. An explicit conversion at a JSON or database boundary is clearer than scattering `.item()` calls through the numerical core.

### Storage and interchange

The `.npy` format saves one array's shape, `dtype`, and data, while `.npz` can bundle several named arrays in one archive. They support exact round trips within NumPy workflows but aren't universal tabular protocols. Choose a format for the consumer before choosing one merely because it is easy to write.

Object arrays may depend on Python pickle. Loading untrusted pickle data can execute code needed to construct objects, so `allow_pickle=True` isn't a general fix for an unfamiliar-file error. Prefer numeric, Boolean, or fixed-string types and keep pickle disabled for untrusted files.

Text formats lose some type and shape information. CSV can't naturally represent an arbitrary-dimensional array or completely preserve floating format, byte order, and structured `dtype`. A text export needs an accompanying contract for columns, missing values, precision, and shape reconstruction.

When exchanging arrays with Pandas, a machine learning framework, or a native extension, check whether conversion copies data, shares writeable memory, or imposes device and byte-order requirements. Similar conversion method names can have different ownership semantics. A reliable interface test changes one side before and after exchange and checks whether the other side should change.

<!-- /deep -->

[Checkpoint: datascience/numpy](https://codewiki.com/datascience/numpy/#checkpoint)

## Further reading

- [NumPy: the absolute basics for beginners](https://numpy.org/doc/stable/user/absolute_beginners.html)
- [NumPy: broadcasting](https://numpy.org/doc/stable/user/basics.broadcasting.html)
- [NumPy: copies and views](https://numpy.org/doc/stable/user/basics.copies.html)
- [NumPy: random generator](https://numpy.org/doc/stable/reference/random/generator.html)
- [NumPy: linear algebra](https://numpy.org/doc/stable/reference/routines.linalg.html)