Which points count?

Which points count?#

Points inside and outside a unit quarter circle

Random points in a square. Count the ones inside the curve.

Code behind the opening picture#

This uses Python lists, random, and math, with Matplotlib for drawing; there are no NumPy array operations. The opening image is embedded so it appears immediately. Run the code below and then display(circle_fig) to redraw it. These are the same first 1,500 seeded points used later.

import random
from math import cos, sin, pi
import matplotlib.pyplot as plt

preview_rng = random.Random(142)
preview_points = [(preview_rng.random(), preview_rng.random())
                  for _ in range(1_500)]
inside_x, inside_y, outside_x, outside_y = [], [], [], []
for x, y in preview_points:
    if x*x + y*y < 1.0:
        inside_x.append(x)
        inside_y.append(y)
    else:
        outside_x.append(x)
        outside_y.append(y)
angles = [i * pi / 400 for i in range(201)]
circle_fig, ax = plt.subplots(figsize=(6.4, 4.3), layout="constrained")
ax.scatter(inside_x, inside_y, s=10, color="#007c91", label="Counted")
ax.scatter(outside_x, outside_y, s=12, marker="x",
           color="#b64f22", label="Not counted")
ax.plot([cos(t) for t in angles], [sin(t) for t in angles], color="#172b4d")
ax.set(xlabel="x", ylabel="y", xlim=(0, 1), ylim=(0, 1), aspect="equal")
ax.legend(loc="lower left")
plt.close(circle_fig)

Make the slow part faster#

Compiled code and Python interop#

Same calculation. Same inputs. Check the answer.

PHYS 142/242 · Lab 08

Presenting this notebook#

Open Reveal Slideshow in DataHub. Use the right arrow to move through the slides. Click a code cell and press Shift+Enter to run it. Run the code in order; each demo uses earlier definitions. Ctrl+Enter runs a selected cell while keeping it selected.

This notebook starts with empty outputs for live teaching. For a rehearsal, restart the kernel and run all cells in ordinary notebook view, then open the slideshow. Clear outputs before saving the handout again. The notebook includes its presentation styling in the RISE metadata.

Cells marked Skip contain reading notes and optional reference material. They remain available in notebook view. The Julia and JuliaCall examples are fenced Markdown, not Python cells, so Run All does not install or launch Julia.

Suggested pacing for the 80-minute lab: Python and NumPy 20 minutes; Numba and timing 20; C++ and Julia 20; student experiments and discussion 20. Adjust to the room.

What we’ll do#

  1. Check a small answer by hand.

  2. Run the same calculation in Python, NumPy, Numba, and C++.

  3. Measure the work and decide what is worth compiling.

We’ll also read the Julia version, then try a change of our own.

Start here#

  • Python 3 on DataHub. Run the code in order.

  • NumPy, Matplotlib, Numba. We’ll import them when we use them.

  • C++ needs g++. Julia is a read-through today.

import sys
import random
from math import sqrt
from statistics import median
from time import perf_counter
from timeit import repeat

print("Python:", sys.version.split()[0])

Setup note#

The Python examples also work in Colab; the C++ example requires Linux with g++. If a Python import is missing, install the needed packages before presenting, then restart the kernel:

%pip install numpy matplotlib numba

If the compiler is missing, ask me before changing your runtime during the walkthrough.

Where does the numerical work happen?#

What we write

What runs the work

Python loop

Python handles each iteration

NumPy expression

Compiled array routines

Numba @njit function

Compiled numerical loop

C++ function called from Python

A library we compile first

Using NumPy already means using compiled code.

One small calculation#

Draw points uniformly in the unit square. Count the ones inside the quarter circle:

\[ x^2 + y^2 < 1, \qquad \widehat{\pi} = 4\,\frac{N_{\mathrm{inside}}}{N}. \]

  • Keep the points fixed when comparing implementations.

  • Exclude the boundary: we chose <, not <=.

Predict the count#

Which of these four points count?

The outer list holds the points. Each tuple holds one (x, y) pair.

small_python = [
    (0.0, 0.0),
    (0.3, 0.4),
    (1.0, 0.0),
    (1.0, 1.0),
]

Lists and tuples#

Lists are mutable: we can append, remove, or replace entries. Tuples are immutable: we cannot replace their entries. Here each tuple is one coordinate pair, and the list collects the pairs. In for x, y in points, Python unpacks each pair into x and y.

Count with Python#

Read the loop, then check your prediction.

def count_python(points):
    count = 0
    for x, y in points:       # unpack one coordinate pair
        if x*x + y*y < 1.0:
            count += 1       # count only points inside
    return count

assert count_python(small_python) == 2
print("Points inside:", count_python(small_python))

Generate the sample once#

Same seed. Same stored points for every implementation.

seed = 142
rng = random.Random(seed)
points_python = []

for _ in range(100_000):
    x = rng.random()
    y = rng.random()
    points_python.append((x, y))

for point in points_python[:3]:
    print(point)

Reading note: list comprehensions#

The explicit loop can also be written as a list comprehension:

rng = random.Random(seed)
points_python = [(rng.random(), rng.random()) for _ in range(100_000)]

Read it as: “for each repetition, make a coordinate pair and put it in the list.” The expression that makes each entry comes first. The _ means we do not need the loop counter; underscores inside 100_000 only make the number easier to read.

Both versions draw two random numbers per point. Reset rng = random.Random(seed) before each version if you want to compare their lists exactly. Use whichever form makes your code easier to read.

Estimate π#

How close should 100,000 random points get us?

count = count_python(points_python)
p_inside = count / len(points_python)
pi_estimate = 4 * p_inside
pi_stderr = 4 * sqrt(p_inside * (1 - p_inside) / len(points_python))

print(f"N = {len(points_python):,}, seed = {seed}")
print(f"pi estimate = {pi_estimate:.5f} +/- {pi_stderr:.5f}")

This is an estimated standard error for independent points.

That uncertainty is an estimated standard error from counting independent hits and misses. It describes the sampling noise in this example. It isn’t an error formula to copy into a correlated simulation.

Two different checks#

  • Physics/statistics check: is the estimate reasonable given its sampling noise?

  • Implementation check: do the functions return the exact same count on the exact same points?

A plausible estimate of π is not enough to prove the implementations agree.

Put the same points in an array#

Convert the coordinates. Don’t draw a new sample.

import numpy as np

points = np.array(points_python, dtype=np.float64)
small = np.array(small_python, dtype=np.float64)

print("Shape:", points.shape, "| dtype:", points.dtype)
print("First point:", points[0])

Choose storage that matches the physics#

Our coordinates are real, so we chose float64. A wavefunction is complex.

Before compiling or passing an array to C++, check what its dtype can represent.

print("Coordinate storage:", points.dtype)
z = np.array([1 + 2j, -3j], dtype=np.complex128)
rotated = z * np.exp(0.4j)

np.testing.assert_allclose(np.abs(rotated), np.abs(z),
                           rtol=1e-12, atol=1e-12)
print("Complex input:", z.dtype, "| Output:", rotated.dtype)

Why check dtype here?#

We just chose the array storage used by the compiled routines. Real coordinates fit in float64; a complex wavefunction needs complex128 (or another appropriate complex dtype). Using only the real part would change the state. A phase rotation is a small independent check: its magnitude should be unchanged.

np.zeros(n) allocates real storage; np.zeros(n, dtype=np.complex128) can store complex values. Later, our C++ interface explicitly requires float64 for the coordinate array. Matching its argument declarations to the data is part of correctness.

assert_allclose allows small rounding differences. The tolerances here are for this tiny check. For the circle example we compare integer counts on fixed inputs, so we use exact equality.

Count with NumPy#

The comparison makes a Boolean array. Count its True entries.

def count_numpy(points):
    radius_squared = points[:, 0]**2 + points[:, 1]**2
    return np.count_nonzero(radius_squared < 1.0)

assert count_numpy(small) == 2
assert count_numpy(points) == count
print("Python and NumPy agree:", count)

Short code; compiled array operations; some temporary arrays.

This is short and readable. It also creates intermediate arrays for the squares, sum, and comparison. That can matter for larger inputs. A compiled loop can do the arithmetic and counting in one pass without keeping those arrays.

That doesn’t tell us which version is faster on this machine. We still need to measure it.

Compile the loop with Numba#

Same arithmetic. Give the function an array and add @njit.

from numba import njit

@njit
def count_numba(points):
    count = 0
    for x, y in points:
        if x*x + y*y < 1.0:
            count += 1
    return count

What belongs inside the compiled function?#

Numba works well with supported numerical loops and arrays. Keep plotting, file handling, and setup outside this function. It cannot compile an arbitrary notebook. Pass changing parameters as arguments rather than relying on edited global variables.

The first call includes compilation#

Predict: will this first call reflect the speed of later calls?

start = perf_counter()
numba_count = count_numba(points)
first_call_seconds = perf_counter() - start

assert numba_count == count
assert count_numba(small) == 2
print(f"First call: {first_call_seconds:.4f} s")
print("Counts agree:", numba_count)

Rerunning this cell reuses the compiled function. A fresh-kernel run measures compilation again.

Numba specializes on argument types. A different array dtype can trigger another compilation; changing only its length usually doesn’t.

For later work, pass changing parameters as arguments. Don’t hide them in a global variable and assume compiled code will notice every edit. Start with ordinary @njit; options such as fastmath=True can change floating-point behavior.

Now compile the same loop ourselves#

Numba compiled Python for us. For C++, we write and build the library.

Python will call it through ctypes. First, check g++.

import ctypes
from pathlib import Path
import subprocess

compiler = subprocess.run(["g++", "--version"], check=True,
                          capture_output=True, text=True)
print(compiler.stdout.splitlines()[0])

1. Write the C++ counter#

Two numbers per row: x is at 2*i; y is at 2*i + 1.

%%writefile circle.cpp
#include <cstddef>
extern "C" std::size_t count_inside(
    const double* points, std::size_t n) {
    std::size_t count = 0;
    for (std::size_t i = 0; i < n; ++i) {
        double x = points[2*i];
        double y = points[2*i + 1];
        if (x*x + y*y < 1.0) {
            ++count;
        }
    }
    return count;
}

Writing the library source#

%%writefile writes this cell to circle.cpp in the current working directory. extern "C" gives the function a name that ctypes can look up. std::size_t is used for both the number of points and the count.

2. Compile, then load#

-O3: optimize. -shared -fPIC: build a loadable library.

subprocess.run([
    "g++", "-O3", "-std=c++17", "-shared", "-fPIC",
    "circle.cpp", "-o", "libcircle.so"
], check=True)

library = ctypes.CDLL(str(Path("libcircle.so").resolve()))

3. Declare the interface#

Check dtype, layout, shape, and the return type.

point_array = np.ctypeslib.ndpointer(
    dtype=np.float64, ndim=2, flags=("C_CONTIGUOUS", "ALIGNED")
)
library.count_inside.argtypes = [point_array, ctypes.c_size_t]
library.count_inside.restype = ctypes.c_size_t

def count_cpp(points):
    if points.ndim != 2 or points.shape[1] != 2:
        raise ValueError("Expected an array with shape (N, 2).")
    return library.count_inside(points, len(points))

Why these checks matter#

The function expects an aligned, C-contiguous float64 array with shape (N, 2). ndpointer checks dtype and memory layout; the wrapper checks the shape. A pointer alone does not carry the array dimensions. The call passes the existing NumPy storage to C++, without copying the whole array. Our C++ function only reads it.

4. Check the C++ answer#

Known small answer → same large count. Then compare all four timings.

assert count_cpp(small) == 2
assert count_cpp(points) == count
print("All four counts agree:", count)

The call passes the existing NumPy storage to C++; it doesn’t copy the whole array. That’s why the layout checks matter. Our function only reads the input.

If you change the C++ source after loading it, restart the kernel and rerun before checking the rebuilt library. Loading the same filename again can reuse the library already in memory.

For your own work, this extra setup is worthwhile when you need an existing C/C++ library or a substantial compiled calculation. For a small numerical Python loop, Numba is often less work to try first.

Compare all four implementations#

Predict the order. Numba is warmed up; C++ is already compiled.

def median_seconds(function, data):
    return median(repeat(lambda: function(data), repeat=5, number=3)) / 3
timings = {
    "Python": median_seconds(count_python, points_python),
    "NumPy": median_seconds(count_numpy, points),
    "Numba": median_seconds(count_numba, points),
    "C++": median_seconds(count_cpp, points),
}
for name, seconds in timings.items():
    speedup = timings["Python"] / seconds
    print(f"{name:8s} {seconds * 1000:9.3f} ms  {speedup:7.1f}x")

What did we actually measure?#

  • Included: counting the same stored points, including function-call overhead.

  • Excluded: imports, sampling, conversion, compilation, and plotting.

  • Different storage: Python uses the list; NumPy, Numba, and C++ use the array.

Does a faster counter make this fixed sample more accurate?

No. It may let us afford more samples. Time the full workflow before claiming the whole program is faster.

Before changing your homework code#

  • Find the slow part. Separate numerical work from plotting.

  • Keep the calculation fixed. Same grid, dtype, algorithm, and sample size.

  • Use existing compiled routines. NumPy products and SciPy solvers already do this.

Keep repeated numerical work inside the function when practical.

Checking an optimization#

Matching implementations catches translation mistakes. Keep an independent physics check and a convergence check too. A slow animation does not imply a slow simulation. Calling across the Python/compiled boundary for every scalar operation can eat up the gain. For correlated simulations, more updates do not automatically give more independent measurements.

Which tool is worth the setup?#

Situation

A useful starting point

Standard array operation or solver

NumPy / SciPy

A custom numerical Python loop

Numba

An existing C/C++ library

A checked Python interface

Check the answer and measure your own workload.

The same loop in Julia#

Read this one. It is not Python syntax.

function count_julia(points)
    count = 0
    for i in axes(points, 1)
        x = points[i, 1]
        y = points[i, 2]
        if x*x + y*y < 1.0
            count += 1
        end
    end
    return count
end

Where Julia fits#

  • Compiled numerical functions, with a first-call cost.

  • Ordinary arrays start at 1. Here the columns are 1 and 2.

  • Same checks still apply. Known answer, same inputs, then timing.

Python can call Julia through JuliaCall. That is an optional setup step for later.

Julia reference: a known answer#

Run this in Julia after defining count_julia:

small = [0.0 0.0; 0.3 0.4; 1.0 0.0; 1.0 1.0]
@assert count_julia(small) == 2

axes(points, 1) gives the row indices. Julia specializes methods for argument types; warm up the function before timing repeated calls. See the earlier Julia notebook for arrays, matrix operations, and compiler output.

Optional reference: JuliaCall#

This optional route has not been verified on our DataHub setup. The attempted install/import did not work during preparation. It is not part of today’s runnable walkthrough.

JuliaCall can call Julia from Python. Installation and the first import can download Julia and set up packages. Prepare and test this in advance if you want to explore it.

Installation command, shown as reference only:

%pip install juliacall

After installation, the Python example is:

from juliacall import Main as jl

jl.seval("sumsq(values) = sum(abs2, values)")
values = np.array([3.0, 4.0])
answer = jl.sumsq(values)
assert answer == 25.0

Check returned values and types, and include conversions in timings when your real workflow needs them. Pass a whole array when practical, rather than crossing the language boundary once per element.

Your turn: change the region#

Give the Python and Numba counters a radius argument.

  1. Predict a count for points inside, outside, and on the boundary.

  2. Replace the condition with x*x + y*y < radius*radius.

  3. Check both implementations on exactly the same inputs.

Does changing the argument after compilation change the result?

# Keep the original functions above for comparison.
# Define count_python_radius(points, radius) and
# count_numba_radius(points, radius), then test them.

Another experiment to try#

  • Change the workload: compare three sample sizes. Report first-call cost separately.

  • Choose a target in your code: what would you time, and what small check would catch a mistake?

Be ready to show the check that supports your result.

You don’t need to rewrite an assignment during lab.

Before you keep an optimization#

What got faster?

What stayed the same?

What check would catch a changed answer?

Quick reference#

If you need to…

Start here

Speed up a custom numerical loop

Put it in a function with explicit inputs; try @njit.

Multiply matrices or solve a standard numerical problem

Check NumPy/SciPy’s existing routines first.

Compare implementations

Use the same inputs; check a tiny known answer before a larger case.

Time compiled code

Separate setup, first call, repeated calls, and total runtime.

Preserve complex values

Check input, output, and preallocated array dtypes.

Call C/C++

Declare argument/return types and check array shape/layout.

Reproduce a stochastic check

Record the seed and reuse the generated inputs when comparing code.

The same seed in different languages doesn’t promise the same random stream. For a direct comparison, pass the same samples. For a larger simulation, compare observables and uncertainties too.

Before you keep an optimization: What got faster? What stayed the same? What check would catch it if you changed the answer?

More detail when you need it#

Contents