Lecture 4 · NumPy, Matplotlib, and Readable Evidence

NumPy, Matplotlib, and readable evidence¶

Arrays, figures, and executable code¶

Use one small synthetic brightness record to inspect arrays and figures.

Run the Python cells in order. They print the arrays and render the figures used in this lecture.

Use the docs for the full API¶

We show a small, useful subset of NumPy and Matplotlib. Use the official docs to check the full API: signatures, options, return values, errors, and examples.

• NumPy: user guide · API reference

• Matplotlib: documentation · API reference · cheatsheets · Figure anatomy

The through-line¶

State what positions, axes, and marks mean before interpreting a result.

Today’s question¶

How do we organize values across positions and axes, select and combine them safely, and make a figure a reader can check?

By the end of class, you should be able to:

• create and inspect NumPy arrays with shape, size, and dtype;

• use indexing, slicing, and Boolean masks while preserving aligned data;

• read and index 2-D and ND arrays by their named axes and shape order;

• use a length-one axis to broadcast an explicitly owned quantity safely;

• write and check a vectorized calculation; and

• construct a readable Matplotlib figure and state what it shows—and what it cannot establish alone.

One synthetic record supplies the values¶

These synthetic brightness measurements let us practice array operations and plots without reproducing a Project 1 analysis. A periodic reference curve appears later to distinguish calculated values from observations.

I · NumPy arrays: quantities, positions, and elementwise operations¶

Question¶

What is an array, and what must be true before one operation can safely apply to many values?

You will leave this section able to¶

• distinguish a Python list from a NumPy array;

• read shape, size, and dtype;

• create a numerical grid with np.linspace; and

• explain elementwise arithmetic.

Through-line focus: named quantity → numerical array.

A list and an array make different promises¶

import numpy as np

brightness_list = [1.000, 0.997, 0.991]
brightness_array = np.array(brightness_list)

print("list * 2:", brightness_list * 2)
print("array * 2:", brightness_array * 2)

list * 2: [1.0, 0.997, 0.991, 1.0, 0.997, 0.991]
array * 2: [2.    1.994 1.982]

Shape says how positions are organized¶

observation_time_days = np.linspace(0.0, 6.0, 13)
observed_brightness = np.array([
    0.986, 0.971, 0.973, 0.984, 0.999, 0.997, 0.985,
    0.973, 0.970, 0.986, 0.997, 0.999, 0.984,
])

print("time:", observation_time_days)
print("brightness:", observed_brightness)
print("time shape:", observation_time_days.shape)
print("brightness shape:", observed_brightness.shape)
print("number of observations:", observed_brightness.size)
print("dtype:", observed_brightness.dtype)

time: [0.  0.5 1.  1.5 2.  2.5 3.  3.5 4.  4.5 5.  5.5 6. ]
brightness: [0.986 0.971 0.973 0.984 0.999 0.997 0.985 0.973 0.97  0.986 0.997 0.999
 0.984]
time shape: (13,)
brightness shape: (13,)
number of observations: 13
dtype: float64

A numerical grid names its endpoints and count¶

time_days = np.linspace(0.0, 4.0, 9)

print("time grid:", time_days)
print("shape:", time_days.shape)
print("spacing:", time_days[1] - time_days[0], "days")

time grid: [0.  0.5 1.  1.5 2.  2.5 3.  3.5 4. ]
shape: (9,)
spacing: 0.5 days

Common array generators¶

integer_steps = np.arange(0, 5)
log_values = np.logspace(0, 3, 4)
blank_by_source = np.zeros((3, 5))
unit_weights = np.ones(5)
missing_values = np.full(5, np.nan)

print("integer steps:", integer_steps)
print("log values:", log_values)
print("blank shape:", blank_by_source.shape)
print("weights:", unit_weights)
print("missing values:", missing_values)

integer steps: [0 1 2 3 4]
log values: [   1.   10.  100. 1000.]
blank shape: (3, 5)
weights: [1. 1. 1. 1. 1.]
missing values: [nan nan nan nan nan]

np.logspace(0, 3, 4) returns 10**0 through 10**3. For a floating-point grid with a stated number of values, use np.linspace rather than np.arange.

Elementwise arithmetic follows positions¶

background_level = 0.003
corrected_brightness = observed_brightness - background_level

print("first three corrected values:", corrected_brightness[:3])
print("shape preserved:", corrected_brightness.shape)

first three corrected values: [0.983 0.968 0.97 ]
shape preserved: (13,)

dtype is a numerical representation, not a unit¶

integer_counts = np.array([1000, 997, 991])
normalized_brightness = integer_counts / integer_counts[0]

print("count dtype:", integer_counts.dtype)
print("normalized dtype:", normalized_brightness.dtype)
print("normalized values:", normalized_brightness)

count dtype: int64
normalized dtype: float64
normalized values: [1.    0.997 0.991]

Reference · array creation and inspection¶

Rule: state what each position represents before combining arrays.

II · Indexing, slicing, and Boolean selection¶

Question¶

How do we retrieve or select values while preserving the relationships among aligned arrays?

You will leave this section able to¶

• retrieve the first, last, or a contiguous range of positions;

• construct a Boolean mask from a stated rule;

• apply the same mask to paired arrays; and

• count selected values before interpreting them.

Through-line focus: paired arrays → same selection.

Indexing starts at 0; negative indices count from the end¶

first_time_days = observation_time_days[0]
third_time_days = observation_time_days[2]
last_time_days = observation_time_days[-1]
penultimate_time_days = observation_time_days[-2]
third_brightness = observed_brightness[2]
early_time_days = observation_time_days[:3]
early_brightness = observed_brightness[:3]

print("first observation:", first_time_days, "days")
print("third observation:", third_time_days, "days;", third_brightness)
print("last observation:", last_time_days, "days")
print("one before last:", penultimate_time_days, "days")
print("early times:", early_time_days)
print("early brightness:", early_brightness)

assert observation_time_days[-1] == observation_time_days[len(observation_time_days) - 1]

first observation: 0.0 days
third observation: 1.0 days; 0.973
last observation: 6.0 days
one before last: 5.5 days
early times: [0.  0.5 1. ]
early brightness: [0.986 0.971 0.973]

A Boolean mask selects the same positions from both arrays¶

low_brightness = observed_brightness < 0.985

print("mask:", low_brightness)
print("selected times:", observation_time_days[low_brightness])
print("selected brightness:", observed_brightness[low_brightness])
print("selected count:", np.count_nonzero(low_brightness))

mask: [False  True  True  True False False False  True  True False False False
  True]
selected times: [0.5 1.  1.5 3.5 4.  6. ]
selected brightness: [0.971 0.973 0.984 0.973 0.97  0.984]
selected count: 6

Arrays need an explicit many-value question¶

try:
    if low_brightness:
        print("This will not run.")
except ValueError as error:
    print(type(error).__name__ + ":", error)

print("any low values?", np.any(low_brightness))
print("all low values?", np.all(low_brightness))

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
any low values? True
all low values? False

A threshold is a selection rule¶

Reference · indexing and Boolean selection¶

Rule: use the same selection on every array that describes the same observations.

III · 2-D and ND arrays: matrices, axes, and shape¶

Question¶

How do we represent values across multiple named axes without losing the meaning of any dimension?

You will leave this section able to¶

• read ndim and a two-dimensional shape;

• state what each axis and position represents;

• index a row, column, or one specific value;

• state what axis=0 and axis=1 leave behind; and

• broadcast one owned value per source across time.

Through-line focus: named positions → named axes.

A 2-D array can represent a labeled matrix¶

source_labels = np.array(["A", "B", "C"])
sample_time_days = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
brightness_by_source = np.array([
    [1.000, 0.998, 0.985, 0.997, 1.001],
    [0.980, 0.979, 0.966, 0.978, 0.981],
    [1.020, 1.018, 1.006, 1.017, 1.021],
])

print("ndim:", brightness_by_source.ndim)
print("shape:", brightness_by_source.shape)
print("row labels:", source_labels)
print("column times (days):", sample_time_days)

ndim: 2
shape: (3, 5)
row labels: ['A' 'B' 'C']
column times (days): [0. 1. 2. 3. 4.]

Matrix order and coordinate order are different conventions¶

For a matrix or ordinary 2-D array, use array[row, column]: first vertical row, then horizontal column.

Cartesian coordinates are conventionally written (x, y): horizontal position first, then vertical position. An image stored in a 2-D array is therefore commonly accessed as image[y, x], not image[x, y].

Shape and indexing use the same axis order¶

source_index = 1
time_index = 2

brightness_by_time_source = brightness_by_source.T

print("axis order:", "(source, time)")
print("shape:", brightness_by_source.shape)
print("value at [source, time]:", brightness_by_source[source_index, time_index])
print("transposed axis order:", "(time, source)")
print("transposed shape:", brightness_by_time_source.shape)
print("same value at [time, source]:", brightness_by_time_source[time_index, source_index])

assert (
    brightness_by_source[source_index, time_index]
    == brightness_by_time_source[time_index, source_index]
)

axis order: (source, time)
shape: (3, 5)
value at [source, time]: 0.966
transposed axis order: (time, source)
transposed shape: (5, 3)
same value at [time, source]: 0.966

Indexing a 2-D array names both axes¶

source_index = 1
time_index = 2

print("one source:", brightness_by_source[source_index])
print("one time across sources:", brightness_by_source[:, time_index])
print(
    "one measurement:",
    source_labels[source_index],
    "at day",
    sample_time_days[time_index],
    "=",
    brightness_by_source[source_index, time_index],
)

one source: [0.98  0.979 0.966 0.978 0.981]
one time across sources: [0.985 0.966 1.006]
one measurement: B at day 2.0 = 0.966

A reduction consumes one axis¶

mean_brightness_per_source = brightness_by_source.mean(axis=1)
mean_brightness_per_time = brightness_by_source.mean(axis=0)

print("per-source shape:", mean_brightness_per_source.shape)
print("per-source means:", mean_brightness_per_source)
print("per-time shape:", mean_brightness_per_time.shape)
print("per-time means:", mean_brightness_per_time)

per-source shape: (3,)
per-source means: [0.9962 0.9768 1.0164]
per-time shape: (5,)
per-time means: [1.         0.99833333 0.98566667 0.99733333 1.001     ]

Broadcasting: one operation across compatible shapes¶

Definition: Broadcasting is NumPy’s shape-matching rule for elementwise operations. NumPy aligns shapes from the right. On each aligned axis, the lengths must match or one length must be 1; missing leading axes behave like length 1. The result uses the larger length on each axis.

NumPy does not need to copy a value into a larger array. Broadcasting describes how the operation is evaluated.

Apply the shape rule before the operation¶

Equal lengths are compatible. A length-one axis can expand. Any other pair of aligned lengths raises ValueError.

See broadcasting on generic rows and columns¶

values_by_row_and_column = np.zeros((3, 5))
per_column = np.array([10, 20, 30, 40, 50])
per_row = np.array([1, 2, 3])

print("per-column result shape:", (values_by_row_and_column + per_column).shape)

try:
    values_by_row_and_column + per_row
except ValueError as error:
    print("per-row without newaxis:", type(error).__name__)

per_row_column = per_row[:, np.newaxis]
print("per-row result shape:", (values_by_row_and_column + per_row_column).shape)
print("per-row result:\n", values_by_row_and_column + per_row_column)

per-column result shape: (3, 5)
per-row without newaxis: ValueError
per-row result shape: (3, 5)
per-row result:
 [[1. 1. 1. 1. 1.]
 [2. 2. 2. 2. 2.]
 [3. 3. 3. 3. 3.]]

Apply broadcasting to one baseline per source¶

baseline_per_source = brightness_by_source[:, 0]

try:
    brightness_by_source / baseline_per_source
except ValueError as error:
    print("without an explicit axis:", type(error).__name__)

baseline_column = baseline_per_source[:, np.newaxis]
normalized_by_source = brightness_by_source / baseline_column

assert normalized_by_source.shape == brightness_by_source.shape
assert np.allclose(normalized_by_source[:, 0], 1.0)
print("baseline shape:", baseline_per_source.shape)
print("explicit baseline shape:", baseline_column.shape)
print("normalized shape:", normalized_by_source.shape)
print("first normalized row:", normalized_by_source[0])

without an explicit axis: ValueError
baseline shape: (3,)
explicit baseline shape: (3, 1)
normalized shape: (3, 5)
first normalized row: [1.    0.998 0.985 0.997 1.001]

Reference · ND arrays, axes, and broadcasting¶

Rule: before combining ND arrays, name each axis, its owner, and the expected output shape.

Practice the 2-D case before using it in larger work¶

The Vectorization and Broadcasting Lab repeats this reasoning with a fresh data table. Complete it after this lecture: predict a shape, write the axis meaning, then run the expression.

IV · Vectorized calculations and shape checks¶

Question¶

How can one stated calculation accept a scalar or many positions while preserving correspondence and checking shape?

You will leave this section able to¶

• write a vectorized function with an explicit input contract;

• check an array result against one scalar result; and

• distinguish a vectorized one-dimensional calculation from an ND-array broadcasting operation.

Through-line focus: paired array → stated calculation.

A vectorized calculation accepts a scalar, list, or array¶

def periodic_reference_brightness(time_days, reference_period_days, dimming_depth):
    """Return a periodic reference curve for finite times in days.

    dimming_depth sets the drop from 1.0 to 1.0 - dimming_depth.
    """
    time_days = np.asarray(time_days, dtype=float)
    if not np.all(np.isfinite(time_days)):
        raise ValueError("time_days must contain only finite values")
    if not np.isfinite(reference_period_days) or reference_period_days <= 0:
        raise ValueError("reference_period_days must be finite and positive")
    if not np.isfinite(dimming_depth) or not 0 <= dimming_depth < 1:
        raise ValueError("dimming_depth must be finite and in [0, 1)")

    phase = 2.0 * np.pi * time_days / reference_period_days
    return 1.0 - dimming_depth * (1.0 + np.sin(phase)) / 2.0


reference_period_days = 3.0
dimming_depth = 0.030
reference_brightness = periodic_reference_brightness(
    observation_time_days,
    reference_period_days,
    dimming_depth,
)

print("reference shape:", reference_brightness.shape)
print("reference values:", np.round(reference_brightness, 4))

reference shape: (13,)
reference values: [0.985 0.972 0.972 0.985 0.998 0.998 0.985 0.972 0.972 0.985 0.998 0.998
 0.985]

Check one vectorized result against one scalar call¶

index = 4
scalar_brightness = periodic_reference_brightness(
    observation_time_days[index],
    reference_period_days,
    dimming_depth,
)

np.testing.assert_allclose(
    reference_brightness[index],
    scalar_brightness,
    rtol=0.0,
    atol=1e-12,
)
print("scalar and vector results agree at index", index)

scalar and vector results agree at index 4

A shape error asks a representation question¶

reference_time_days = np.linspace(0.0, 6.0, 400)

try:
    observation_time_days - reference_time_days
except ValueError as error:
    print(type(error).__name__ + ":", error)

ValueError: operands could not be broadcast together with shapes (13,) (400,) 

Reference · vectorization and compatible shapes¶

Rule: write the shape and physical meaning of every position before combining arrays.

V · Matplotlib: readable figures and bounded evidence¶

Question¶

What must a reader identify to interpret a figure and assess a bounded claim?

You will leave this section able to¶

• create a figure and axes object;

• distinguish observations, a reference calculation, and a threshold;

• choose scales and axis direction deliberately; and

• write a bounded evidence statement from a figure.

Through-line focus: stated calculation → readable figure.

Run this: render the observations¶

fig is the whole figure. ax is the plotting area that owns the data, labels, limits, legend, and grid.

Use fig, ax = plt.subplots() as the default course pattern. It keeps later choices attached to a named axes object and extends cleanly to multi-panel figures.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(
    figsize=(8.5, 4.8),
    layout="constrained",
)
ax.scatter(
    observation_time_days,
    observed_brightness,
    color="#0072B2",
    marker="o",
    label="synthetic observations",
)
ax.set(
    xlabel="Observation time (days)",
    ylabel="Normalized brightness (dimensionless)",
    title="Synthetic observations",
)
ax.grid(alpha=0.25)
ax.legend()
fig.savefig(
    "lecture-04-synthetic-observations.png",
    dpi=200,
    bbox_inches="tight",
)
plt.show()

Read the rendered output: each circle is one observation. No curve is shown yet.

Add a separately named reference curve¶

A line can look more authoritative than the data. Name it as a reference calculation and use a denser, separately named grid so the reader can distinguish calculated values from observations.

reference_time_days = np.linspace(0.0, 6.0, 400)
smooth_reference_brightness = periodic_reference_brightness(
    reference_time_days,
    reference_period_days,
    dimming_depth,
)

comparison_fig, comparison_ax = plt.subplots(
    figsize=(8.5, 4.8),
    layout="constrained",
)
comparison_ax.scatter(
    observation_time_days,
    observed_brightness,
    color="#0072B2",
    marker="o",
    label="synthetic observations",
    zorder=3,
)
comparison_ax.plot(
    reference_time_days,
    smooth_reference_brightness,
    color="#D55E00",
    linestyle="-",
    label="periodic reference curve",
)
comparison_ax.set(
    xlabel="Observation time (days)",
    ylabel="Normalized brightness (dimensionless)",
    title="Synthetic observations and a reference calculation",
    xlim=(0.0, 6.0),
)
comparison_ax.grid(alpha=0.25)
comparison_ax.legend()
comparison_fig.savefig(
    "lecture-04-observations-reference.png",
    dpi=200,
    bbox_inches="tight",
)
plt.show()

Read the rendered output: markers are 13 observations; the solid line uses a denser reference grid.

Add a threshold and complete the figure¶

import matplotlib.pyplot as plt

reference_time_days = np.linspace(0.0, 6.0, 400)
smooth_reference_brightness = periodic_reference_brightness(
    reference_time_days,
    reference_period_days,
    dimming_depth,
)
threshold_brightness = 0.985

brightness_fig, ax = plt.subplots(figsize=(8.5, 4.8), layout="constrained")
ax.scatter(
    observation_time_days,
    observed_brightness,
    color="#0072B2",
    marker="o",
    label="synthetic observations",
    zorder=3,
)
ax.plot(
    reference_time_days,
    smooth_reference_brightness,
    color="#D55E00",
    linestyle="-",
    label="periodic reference curve",
)
ax.axhline(
    threshold_brightness,
    color="#4D4D4D",
    linestyle="--",
    label="illustrative threshold",
)
ax.set(
    xlabel="Observation time (days)",
    ylabel="Normalized brightness (dimensionless)",
    title="Synthetic brightness observations and a reference calculation",
    xlim=(0.0, 6.0),
)
ax.grid(alpha=0.25)
ax.legend()
brightness_fig.savefig(
    "lecture-04-synthetic-brightness-figure.png",
    dpi=200,
    bbox_inches="tight",
)
plt.show()

• Circles are the synthetic observations.

• The solid line is a reference calculation, not a fit.

• The dashed line is an illustrative threshold.

• Marker shape and line style work even without color.

Read the rendered output: markers, the solid line, and the dashed line represent different things. The legend names them.

Use figure anatomy as a checklist¶

from pathlib import Path

export_path = Path("lecture-04-synthetic-brightness-figure.png")
brightness_fig.savefig(
    export_path,
    dpi=200,
    bbox_inches="tight",
)
print(f"saved: {export_path} (exists: {export_path.exists()})")

saved: lecture-04-synthetic-brightness-figure.png (exists: True)

Check the axes before exporting¶

print("x label:", ax.get_xlabel())
print("y label:", ax.get_ylabel())
print("title:", ax.get_title())
print("legend labels:", ax.get_legend_handles_labels()[1])

x label: Observation time (days)
y label: Normalized brightness (dimensionless)
title: Synthetic brightness observations and a reference calculation
legend labels: ['synthetic observations', 'periodic reference curve', 'illustrative threshold']

The figure object exposes labels and legend entries for code-based checks.

Scale changes what a reader can see¶

distance_pc = np.array([1.0, 10.0, 100.0, 1_000.0])
relative_signal = distance_pc ** -2

scale_fig, axes = plt.subplots(1, 2, figsize=(10.0, 3.8), layout="constrained")

for ax, title in zip(axes, ("Linear axes", "Log--log axes")):
    ax.plot(distance_pc, relative_signal, marker="o", linestyle="-")
    ax.set(
        xlabel="Distance (pc)",
        ylabel="Relative signal (dimensionless)",
        title=title,
    )
    ax.grid(alpha=0.25)

axes[1].set_xscale("log")
axes[1].set_yscale("log")
scale_fig.savefig(
    "lecture-04-scale-comparison.png",
    dpi=200,
    bbox_inches="tight",
)
plt.show()

The values are positive, so logarithmic axes are allowed. Equal ratios occupy equal distances on a log axis. A log scale changes the view, not the values.

Axis direction is a stated convention¶

lookback_time_gyr = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
normalized_quantity = np.array([1.00, 0.93, 0.87, 0.80, 0.74])
detector_row = np.array([1, 2, 3, 4, 5])
row_signal = np.array([12, 9, 7, 10, 11])

direction_fig, axes = plt.subplots(1, 2, figsize=(10.0, 3.8), layout="constrained")

axes[0].plot(lookback_time_gyr, normalized_quantity, marker="o", linestyle="-")
axes[0].invert_xaxis()
axes[0].set(
    xlabel="Lookback time (Gyr; present at right)",
    ylabel="Normalized quantity (dimensionless)",
    title="An inverted x-axis",
)
axes[0].grid(alpha=0.25)

axes[1].plot(row_signal, detector_row, marker="s", linestyle="-")
axes[1].invert_yaxis()
axes[1].set(
    xlabel="Detector signal (counts)",
    ylabel="Detector row (1 at top)",
    title="An inverted y-axis",
)
axes[1].grid(alpha=0.25)
direction_fig.savefig(
    "lecture-04-axis-direction.png",
    dpi=200,
    bbox_inches="tight",
)
plt.show()

Call invert_xaxis() or invert_yaxis() only when the label states the convention.

State what the figure shows¶

Use three parts:

1. Visible: Several synthetic observations and portions of the reference curve lie below the 0.985 illustrative threshold.

2. Evidence: The labeled axes, circular markers, solid line, and dashed threshold distinguish the three visual elements.

3. Open: The reference calculation has no measured uncertainty, calibration, or competing physical explanation.

Reference · a readable Matplotlib figure¶

Rule: name variables, units, calculation status, and selection rules.

Documentation for the next question¶

Exit check¶

A teammate writes:

ax.plot(observation_time_days, observed_brightness)
ax.axhline(0.985)

Name three additions needed before this can serve as a course figure, and say why each matters. Then answer:

1. Why must one Boolean mask select both observation_time_days and observed_brightness?

2. What does a threshold line show directly? What does it not establish?

3. What does reference_brightness.shape == observation_time_days.shape check, and what does it not prove about the pairing?

4. brightness_by_source has shape (3, 5) and baseline_per_source has shape (3,). Why does their division fail? What does baseline_per_source[:, np.newaxis] change, and why does (3, 1) work?

5. Why might a Cartesian point (x, y) be accessed as image[y, x] in a 2-D image array?

Carry forward¶

state quantities and units

→ build and inspect arrays

→ select with a documented rule

→ evaluate a stated calculation

→ make a labeled figure

→ claim only what the figure supports

Reuse this workflow in later work. The next scientific question determines the data, assumptions, and validation it needs.