03 · Arrays are the working language¶

Scientific Python becomes useful when you stop treating an array as “many numbers” and start treating its shape as part of the model.

Decision: does each axis retain a declared physical meaning?

Output: a vectorised, shape-preserving numerical transformation.

Frame¶

A time series may have shape (500,). A field sampled on a grid might have shape (256, 128). A parameter sweep can have shape (n_oh, n_bo, n_time).

Those axes are not interchangeable.

energy.shape
time.shape
field.shape

Check shapes before values when an operation surprises you.

Predict¶

For

time = np.linspace(0.0, 1.0, 6)
height = time**2

predict:

each shape;
the first and last value;
the exact derivative;
which samples have neighbours on both sides.

Then run it.

Indexing answers a question¶

late = time > 0.5
late_energy = energy[late]

The Boolean array late has the same shape as time. It represents a predicate for every sample. Applying it to a different-length array is a model error, not merely an indexing inconvenience.

Slices preserve order:

interior = height[1:-1]
left = height[:-2]
right = height[2:]

These three arrays have the same shape. That alignment is the basis of a centred difference:

derivative[1:-1] = (right - left) / (2 * step)

Broadcasting needs an axis story¶

oh = np.array([0.001, 0.01, 0.1])
bo = np.array([0.001, 0.01, 0.1, 1.0])

oh * bo fails because shapes (3,) and (4,) do not align. If you want every pair:

oh_grid = oh[:, None]   # shape (3, 1)
bo_grid = bo[None, :]   # shape (1, 4)
combined = oh_grid * bo_grid  # shape (3, 4)

The added axes express the intended outer product. Do not add None until you can say what each axis means.

Vectorisation is not a moral virtue¶

Vectorisation is excellent when:

one operation applies independently to an array region;
shapes express the algorithm;
NumPy performs the work more clearly and efficiently.

A loop is better when:

each case may fail differently;
processing has side effects;
the next step depends on a convergence history;
clarity would be lost in a dense expression.

In CoMPhy post-processing, vectorise arithmetic within one case and loop or parallelise across independent cases. That boundary often matches the scientific structure.

Numerical arrays have a data type¶

values.dtype

Integer division produces floats, but storing into an integer array can still truncate. Strings loaded from a CSV do not become physical numbers until you convert and validate them. nan and inf are floats, so type alone does not guarantee validity.

if not np.all(np.isfinite(values)):
    raise ValueError("values must be finite")

A visible finite-difference contract¶

The reference implementation deliberately accepts only a uniform grid:

"""Numerical reductions with explicit assumptions."""

from typing import Any

import numpy as np
from numpy.typing import ArrayLike, NDArray

from .io import SimulationLog


def _one_dimensional(
    name: str,
    values: ArrayLike,
) -> NDArray[np.float64]:
    array = np.asarray(values, dtype=float)
    if array.ndim != 1:
        raise ValueError(f"{name} must be one-dimensional")
    if not np.all(np.isfinite(array)):
        raise ValueError(
            f"{name} must contain only finite values",
        )
    return array


def central_difference(
    x: ArrayLike,
    y: ArrayLike,
) -> NDArray[np.float64]:
    """Differentiate sampled data on a uniformly spaced grid.

    Interior points use a centred difference. The endpoints use a
    second-order one-sided difference. A non-uniform grid is
    refused so the numerical assumption remains explicit.
    """

    coordinate = _one_dimensional("x", x)
    values = _one_dimensional("y", y)
    if coordinate.size != values.size:
        raise ValueError("x and y must have the same length")
    if coordinate.size < 3:
        raise ValueError("at least three samples are required")

    spacing = np.diff(coordinate)
    if np.any(spacing <= 0):
        raise ValueError("x must be strictly increasing")
    if not np.allclose(
        spacing,
        spacing[0],
        rtol=1e-10,
        atol=1e-14,
    ):
        raise ValueError(
            "central_difference requires a uniform grid",
        )

Refusing non-uniform spacing is better than applying a uniform formula silently. A more general algorithm can be added when the course needs it.

Verify¶

Complete Exercise 03. Use \(y=x^2\) because the second-order boundary and interior formulas should reproduce its derivative to roundoff on a uniform grid.

Then change one coordinate to make the grid non-uniform. The correct behaviour for this exercise is a clear error.

Reflect¶

For every array in your solution, annotate:

name:
shape:
axis meaning:
units:
valid range:

If that feels excessive for a five-line exercise, imagine the same ambiguity inside a three-dimensional parameter sweep.