C++ for Data Scientists: Extend Python with pybind11

C++ belongs in a data-science workflow when a measured numerical kernel has become both important and stable. It is not the natural first language for loading a CSV, choosing a model or inspecting residuals. Python and R already make those tasks legible. Native code becomes defensible when profiling isolates a custom calculation, the calculation has a precise statistical contract, and its callers need an extension rather than a separate service. This project is an applied extension in Programming for Data Scientists.

This guide moves one such calculation behind a Python boundary. It consumes the row-level predictions.csv produced by the Python modelling guide and checks its result against the independent r_audit.csv from the R audit. Python still owns file parsing, schema validation and the random resampling plan. C++ computes a paired day-block bootstrap from already validated NumPy arrays. A readable Python oracle, negative boundary tests and sanitizers determine whether the extension is trustworthy.

Start with a native-extension decision, not C++ syntax

Before creating a binding, profile the real workflow with representative input. The decision is not “could this loop be written in C++?” Nearly every loop could be. Ask whether the measured cost is concentrated in code that NumPy, SciPy, Numba or an existing native library does not already execute efficiently. Then ask whether the kernel’s inputs, outputs and numerical meaning are stable enough to support a compiled interface.

A native extension adds obligations: a compiler matrix, Python ABI compatibility, wheel or local-build policy, memory-safety review, sanitizer coverage and a second debugging stack. Those costs are worthwhile only when the boundary stays small. For this learning project the kernel is intentionally narrow:

paired_day_bootstrap(
  actual, baseline_prediction, model_prediction,
  day_offsets, sample_plan
) -> (point_estimate, replicates)

The project makes no speed claim from the 2,168-row historical test set. It is a teaching workload for interface discipline, not a performance benchmark. If a vectorised Python oracle already meets the operational requirement, keeping that implementation is usually the better engineering choice.

Preserve the analytical handover

The extension does not load the original Bike Sharing file and does not retrain the ridge model. The Python stage has already fixed the chronological holdout and exported one record per observed hour. The R stage has independently recalculated the errors and used calendar days as bootstrap blocks.

For each row, the paired absolute-error difference is

d_i = abs(actual_i - model_prediction_i)
      - abs(actual_i - baseline_prediction_i)

Negative values favour the candidate model. Pairing matters because both predictions are evaluated against the same observation; resampling their errors independently would discard that covariance. Whole-day blocks matter because adjacent hourly observations are not plausibly independent. A replicate samples test days with replacement, concatenates every difference in the selected blocks and averages the resulting values.

r_audit.csv is evidence, not an input that tells the native function what answer to return. The wrapper checks that the native point estimate agrees with the independently stored R estimate within the declared tolerance. It does not expect C++ and R bootstrap bounds to be byte-identical: their random-number generators and sample plans differ. Instead, Python supplies one deterministic plan to both the readable oracle and the native kernel, so those two implementations can be compared replicate by replicate.

Use a published C++ baseline

This package targets C++23, the published language baseline available during this guide’s review. C++26 remained a working draft on 16 August 2026; draft features are therefore excluded from the reproducibility contract. That distinction prevents an evergreen title from implying that every compiler already implements an unpublished standard.

The build pins GCC 14.2.0 on the stable Linux runner. A compiler version is only one part of the environment: Python 3.14.7, CMake 4.4.2, Ninja 1.13.2, clang-format 18.1.8, pybind11 3.0.4 and NumPy 2.5.2 are explicit too. The Python-distributed wheels are hash-locked. Because the ninja Python distribution does not publish version 1.13.2, a separate toolchain lock validates the official release asset’s URL, byte count, SHA-256 digest and archive member before installation. Configure only after those checks:

cd examples/programming-for-data-science/cpp
python -m pip install --require-hashes --requirement requirements.lock.txt
python tools/install_ninja.py --destination build/tools
export PATH="$PWD/build/tools:$PATH"
clang-format --dry-run --Werror src/paired_bootstrap.cpp
PYTHONDONTWRITEBYTECODE=1 python -m compileall -q python tests tools
CXX=g++-14 cmake --preset release
cmake --build --preset release
ctest --preset release --output-on-failure

An editor may use clangd or another language server, but editor metadata is not the build definition. CMake and its presets are. This separation lets a reviewer reproduce the package without inheriting an author’s IDE state.

Build the Python module with modern CMake

The key configuration declares the language level and asks CMake for the exact Python interpreter and extension-module headers. NumPy remains pinned in the Python environment; pybind11 then creates a Python-loadable module:

cmake_minimum_required(VERSION 4.4.2)
project(cpp_bootstrap_audit LANGUAGES CXX)

find_package(Python 3.14.7 EXACT
  COMPONENTS Interpreter Development.Module REQUIRED)
set(PYBIND11_FINDPYTHON ON)
find_package(pybind11 3.0.4 EXACT CONFIG REQUIRED)

pybind11_add_module(_paired_bootstrap MODULE src/paired_bootstrap.cpp)
target_compile_features(_paired_bootstrap PRIVATE cxx_std_23)

find_package(Python ...) is preferable to manually assembling include and library paths. It keeps the interpreter and development artefacts aligned. The preset also chooses Ninja, emits a separate build tree and enables strict diagnostics. The native package remains independent of the Astro website build; publishing a programming guide must not make every site contributor install a C++ compiler.

The first useful “hello world” is therefore not console output. It is an import smoke test that proves the module was built for the active interpreter:

PYTHONPATH="build/release:python" python -c \
  "from _paired_bootstrap import paired_day_bootstrap; print(paired_day_bootstrap.__name__)"

An import failure should stop the workflow before any analytical artefact is written.

Let Python own tabular validation

CSV parsing and domain checks remain in Python because that code is easier to inspect beside the existing pipeline. The wrapper requires the exact prediction schema, strictly increasing timestamps, unique instant values, the fixed test window, finite observations and predictions, and contiguous rows for each date. It then creates arrays with explicit native types:

actual_array = np.array(actual, dtype=np.float64, order="C")
baseline_array = np.array(baseline, dtype=np.float64, order="C")
model_array = np.array(model, dtype=np.float64, order="C")

day_offsets = np.array(offsets, dtype=np.int64, order="C")
rng = np.random.Generator(np.random.PCG64(20260815))
sample_plan = rng.integers(
    0,
    len(day_offsets) - 1,
    size=(2000, len(day_offsets) - 1),
    dtype=np.int64,
)

day_offsets is a half-open index: it begins at zero, ends at the row count and increases strictly. Block j spans day_offsets[j]:day_offsets[j + 1]. The two-dimensional sample_plan contains only valid block indices. Generating it outside C++ makes the randomness visible, deterministic and testable without binding a particular C++ random-number engine into the statistical contract.

Conversion must be explicit. np.array(..., order="C") creates named, owned buffers that stay alive across the native call. The wrapper never hands C++ a temporary view with an unclear lifetime. The extension also refuses object arrays, implicit float32 inputs and non-finite values rather than quietly coercing them. Its public contract forbids concurrent input mutation during validation and snapshotting; this package satisfies that rule with fresh private arrays. Once the explicit native snapshot is complete, the GIL-free calculation no longer reads caller-owned buffers.

Make the pybind11 boundary fail closed

The binding accepts exact-dtype, naturally aligned, C-contiguous arrays and performs its own defensive validation. Python-side checks improve error messages, but native code cannot assume every caller used the public wrapper. A shortened signature is:

py::tuple paired_day_bootstrap(
    const py::array& actual,
    const py::array& baseline_prediction,
    const py::array& model_prediction,
    const py::array& day_offsets,
    const py::array& sample_plan);

The implementation rejects arrays unless the three row vectors are one-dimensional and equally sized, offsets are one-dimensional, and the plan is two-dimensional. It checks every floating-point input with std::isfinite, verifies offset endpoints and monotonicity, and bounds-checks every sampled day. It also rejects non-finite paired differences or aggregates, because finite operands near the representable limit can still overflow during subtraction or summation. Empty blocks, zero replicates and integer indices outside the block range are errors.

Those checks are part of the public interface. pybind11’s forcecast flag is deliberately absent because accepting an arbitrary dtype through an invisible copy would hide a producer change. A transposed, strided or byte-offset unaligned buffer also fails instead of being silently materialised or dereferenced through an invalid native pointer. Callers that want conversion must do it visibly in Python.

Use spans and RAII to make lifetimes reviewable

After validation, the binding copies the checked input buffers into native-owned std::vector snapshots. The kernel uses non-owning std::span views over those snapshots, while the replicate output is owned by a pybind11 NumPy array returned to Python. There is no manual new or delete, and no raw pointer escapes its buffer lifetime.

const auto actual_view = std::span<const double>{
    actual_snapshot.data(), actual_snapshot.size()
};

py::array_t<double> output(sample_plan.shape(0));
auto output_view = std::span<double>{
    output.mutable_data(), static_cast<std::size_t>(output.size())
};

RAII is useful because array ownership, temporary vectors and exception cleanup remain deterministic. Smart pointers are not automatically superior: this kernel needs value types and borrowed spans, not heap-allocated object graphs. Object-oriented hierarchies would add indirection without expressing the calculation more clearly. A small pure function with a checked boundary is easier to compare with the Python oracle.

Release the GIL only around pure computation

Validation that reads pybind11 objects happens while the Python Global Interpreter Lock is held. The binding then snapshots every checked input into native-owned storage and allocates the output before releasing the lock around the paired-difference, day-summary, point and replicate calculations:

validate_inputs(actual, baseline, model, day_offsets, sample_plan);
auto snapshots = snapshot_validated_inputs(
    actual, baseline, model, day_offsets, sample_plan);
auto output = allocate_owned_output();
{
    py::gil_scoped_release release;
    compute_differences_day_summaries_point_and_replicates(
        snapshots, output);
}
return py::make_tuple(point_estimate, std::move(output));

Code inside that scope must not create Python objects, raise Python exceptions, resize NumPy arrays or touch reference counts. Native exceptions are captured only after the GIL is reacquired. Releasing the GIL enables other Python threads to run; it is not itself a performance claim and does not justify adding internal parallelism.

The first version keeps the kernel serial. Parallel execution can change floating-point reduction order, sanitizer behaviour and failure reporting. It should be considered only after measurement on the intended workload and then introduced with a stated determinism policy. “C++ supports concurrency” is not a reason to make a statistical audit concurrent.

Prove equivalence with a Python oracle

The oracle is intentionally straightforward. It calculates the row-wise paired difference once, applies the same integer plan, concatenates the selected day slices and takes each replicate mean. The test compares the point estimate and every replicate with a tight absolute tolerance and zero relative tolerance:

native_point, native_replicates = paired_day_bootstrap(
    actual, baseline, model, day_offsets, sample_plan
)
oracle_point, oracle_replicates = python_oracle(
    actual, baseline, model, day_offsets, sample_plan
)

np.testing.assert_allclose(native_point, oracle_point, rtol=0.0, atol=1e-12)
np.testing.assert_allclose(
    native_replicates, oracle_replicates, rtol=0.0, atol=1e-12
)

Tests include the real workload’s uneven block lengths, repeated sampled days and an exactly calculated constant-error case. Negative cases pass float32, non-contiguous and non-finite data; mismatched row counts; invalid offsets; a rank-one plan; and sample indices outside the day range. A good extension test suite spends at least as much attention on rejected inputs as on the happy path because boundary ambiguity is where native code is most dangerous.

Run sanitizers before trusting numerical agreement

Numerical equivalence cannot detect every memory error. The sanitizer preset builds a separate debug configuration with AddressSanitizer and UndefinedBehaviorSanitizer, then executes the C++ and Python-facing tests through CTest:

cd examples/programming-for-data-science/cpp
CXX=g++-14 cmake --preset sanitizers
cmake --build --preset sanitizers
ctest --preset sanitizers --output-on-failure

Compiler warnings are treated as errors. The release preset then rebuilds without borrowing objects from the sanitizer tree and runs the same contract. This is more useful than a generic “compiled successfully” check: it exercises real arrays, exceptions, output ownership and the exact module import path.

Cross-platform wheels are outside this first scope. The stable review target is one documented Linux toolchain. Adding macOS, Windows or multiple Python minor versions would require a deliberate support matrix and artefact policy rather than an untested promise in prose.

Publish deterministic evidence

The command-line wrapper reads predictions.csv and r_audit.csv, creates the deterministic plan, runs both implementations and writes cpp_bootstrap_audit.json. The report records input digests, test-row and day-block counts, the exact toolchain, replicate count, fixed seed, native estimate and percentile bounds, maximum oracle difference, tolerances and agreement with the R point estimate.

JSON keys have a fixed order, floating-point values use one declared serialisation policy, and the file ends with a newline. CI saves the first output, rewrites the report through the same atomic path and compares the two byte sequences. Determinism is evidence that the same reviewed inputs and plan produced the same report; it is not evidence that the historical estimate will generalise.

The underlying UCI sample covers Washington, DC rentals in 2011–2012. It lacks station capacity, bicycle availability, rebalancing operations and prices. The day-block interval describes predictive error differences in the fixed 2012 holdout; it is not current transport evidence, a causal effect, a deployment guarantee or a learning benchmark for language speed.

Decide whether to keep the extension

Keep this C++ boundary only if it earns its maintenance cost. The useful evidence is a production-representative profile, a stable statistical definition, a small typed interface, complete oracle agreement, clean sanitizers and a build matrix that the team can sustain. Remove it if the kernel changes frequently, a maintained library already solves the problem, the data-copy cost dominates, or the Python implementation satisfies the operational requirement.

That decision is the central C++ skill for a data scientist. Modern syntax, modules, templates and parallel algorithms are valuable when a problem needs them, but they do not replace analytical provenance or validation. Begin with the model handover, isolate one measured kernel, make every conversion explicit and require native code to prove that it preserves the result.

Primary references

• The C++ Standards Committee publishes working-group material; the C++26 editors’ report identifies the draft used during this review.
• The pybind11 documentation covers NumPy arrays and its CMake helpers.
• The CMake 4.4 release notes, Python 3.14.7 release and GCC 14 changes document the review toolchain.
• The UCI Bike Sharing record provides the source description, DOI and licence for the historical teaching dataset.